{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths for loading/storing results.\n",
    "name = 'results/Gaussian_MINE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.gaussian import Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 1000   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "\n",
    "rep = 1             # number of repeated runs\n",
    "d = 6               # number of dimensions for X (and Y)\n",
    "\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "g = Gaussian(sample_size=sample_size,rho=rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = g.data\n",
    "        X[i,:,j] = data[:,0]\n",
    "        Y[i,:,j] = data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df7SdZXXnvztcco2BVG9hGlHujY2U1sm6S20GZY1VplJDCtWW1o6ty9hx1Szq2KldmjglKCqhC8gsll214zIzUon1FzPE6iLQoONQxmmwJi6MQXCGWG5ACKDXNBAzide7549z9r37PPd5f53znvO+57zfz1pZ3HPP+77P8wbY+3n23s93i6qCEEJI81hW9QQIIYRUAx0AIYQ0FDoAQghpKHQAhBDSUOgACCGkodABEEJIQ6EDIANHRO4RkT8c0Fh/JCJPisizIvKzOa5/REQuHcTc6oCIfFJEtlc9D1INdACkL7QN6cm24X1SRP5aRM4q+Iw1IqIiMtblHM4EcDOA16vqWar6w26ek/J8FZGXlPlMQgYJHQDpJ7+hqmcBeAWAfwXgmgGP/3MAngPggQGPS8hQQAdA+o6qfh/AXQDWhd+JyDIRuUZEZkTkKRHZJSI/0/763vY/j7V3EhdH7h8XkY+IyOPtPx9p/+4XAHzX3f/V2NxE5K3tsX8oItuC7y4SkX0ickxEnhCRj4rI8vZ3Nrdvtef2b0Xk+SJyh4g8LSI/av/8oqS/FxF5n4h8X0SeEZHvisjrssZtf68i8k4R+b/te68TkbXte46LyG1unpeIyGMicrWI/KC9M3tLypyuEJH722P/g4hMZ82XDDGqyj/8U/ofAI8AuLT98/lorcKva3++B8Aftn9+O4CHAfw8gLMA7AbwqfZ3awAogLGUcT4M4D4A/wLAuQD+wY2Tej+AlwJ4FsBrAIyjFS6ac/P+ZQCvAjDWftaDAN7t7lcAL3GffxbAbwN4LoCzAfw3AH+bMPaFAB4FcJ6b69oC434JwCoA/xLAKQD/o/13+DMAvgPgbe1rL2m/083td3wtgBMALmx//0kA29s/vwLAUwBeCeAMAG9r/3scT5sv/wzvH+4ASD/5WxE5BuBrAP4ewJ9HrnkLgJtV9Xuq+iyAPwPw5gJx/7cA+LCqPqWqTwP4EIC35rz3dwDcoar3quopAO8HMG9fquoBVb1PVedU9REAH0fLgEZR1R+q6u2q+mNVfQbA9SnX/xQtw/pSETlTVR9R1cMFxr1RVY+r6gMADgG4u/13+M9o7bZeHlz/flU9pap/D2APgN+NzOkdAD6uql9X1Z+q6q1oOZdXpc2XDC90AKSf/KaqPk9Vp1T1nap6MnLNeQBm3OcZtFa+P5dzjNj95xW491H7oKonACwkikXkF9phnKMichwtB3ZO0sNE5Lki8vF2SOk4WiGs54nIGeG1qvowgHcD+CCAp0TkcyJyXoFxn3Q/n4x89gn3H7XfzUj6O5oC8J52+OdY23mfj9aqP3G+ZHihAyBV8zhahseYRCtk8SRaoY5u7n8859hPoGXgALQMOFphHONjAB4CcIGqrgJwNQBJed570AqVvLJ9/Wvs0bGLVfUzqvrq9vwVwI1djpvF80Vkpfuc9Hf0KIDr207b/jxXVT+bMV8ypNABkKr5LIA/FZEXt8tE/xzA51V1DsDTaIVkfj7j/mtE5FwROQfABwD8Tc6x/zuAK0Tk1e2k6YfR+f/E2QCOA3hWRH4RwB8F9z8ZzO1stFbfx0RkAsC1SQOLyIUi8qsiMg7g/7Xv+2nOcbvhQyKyXER+BcAVaOUnQv4LgKtE5JXSYqWIXC4iZ2fMlwwpdACkam4B8Cm0wiX/hJZx+WMAUNUfoxVH/9/tkMSrIvdvB7AfwEEA3wbwzfbvMmnHz/89gM+gtRv4EYDH3CXvBfD7AJ5Byzh+PnjEBwHc2p7b7wL4CIAVAH6AVmL671KGHwdwQ/vao2glsa/OOW5RjqL1bo8D+DSAq1T1ofAiVd2PVh7go+3rHwbwBznmS4YUUWVDGEJGFRG5BMDfqGpiOSppLtwBEEJIQ6EDIISQhsIQECGENBTuAAghpKF0pbJYFeecc46uWbOm6mkQQshQceDAgR+o6rnh74fKAaxZswb79++vehqEEDJUiMhM7PcMARFCSEOhAyCEkIZCB0AIIQ2FDoAQQhoKHQAhhDQUOgBCCGkodACEEDIAtuzahy279lU9jQ7oAAghpKEM1UEwQggZNmzVf3BmtuPzjk0XVzYno7IdgIg8R0T+UUS+JSIPiMiHqpoLIYQ0kSp3AKcA/KqqPisiZwL4mojcpar3VTgnQggpFVvp12nlb1TmALSlQ/1s++OZ7T/UpiaE1J46GvNuqDQHICJnADgA4CUA/kpVvx65ZjOAzQAwOTk52AkSQkhJ1NFZ1KIhjIg8D8AXAPyxqh5Kum79+vVKNVBCSFWECd3pqQkA9TTuHhE5oKrrw9/XogxUVY8BuAfAZRVPhRBCGkNlISARORfAT1T1mIisAHApgBurmg8hhGRR54RuN1SZA3gBgFvbeYBlAG5T1TsqnA8hhDSKKquADgJ4eVXjE0JItwz7yt+oRQ6AEEKGiTrq+nQDHQAhhDQUagERQkhO6qzr0w3cARBCGkeREM6ohHticAdACCE5YRkoIYQMKUVCOKMW7olBB0AIIQUZFSdQCy2gvFALiBDSK1t27cPho8exdvWqXIZ8FFb+tdYCIoQQEqefSWiGgAghlTHI1XUY07ffZY09zCv/LOgACCGkhgwiCU0HQAgZON0YN39NN8Zw1Eo4y4AOgBBCasggHBYdACFk4BQxbuFu4cqb9uLEqbnc9yeNnUUTdgp0AISQxjIMRr6fc6MDIIRURh7jFtstdGu4u9lx9MNJ1MXx0AEQQhpHE2Qe8sCTwISQkSc0+CvHW2tfyyVMT00A6HQA/Vz52zxi4/aDpJPA3AEQQhrH2tWrOj43beVvcAdACGkM4aq+qtDPoMflDoAQ0heG2Yg2deVv0AEQQkqjbsnUcD51mVdd5kEHQAhJpGjDFJNZHsSc/MEwANi9dUNfxx1FKnMAInI+gF0AVgOYB7BTVf+iqvkQQnrj8NHjOHFqDgdnZivfCbDMMx9V7gDmALxHVb8pImcDOCAiX1bV71Q4J0II8hlQn0i1lb+XWu4XNq6t/HuRhYiR9ZxRciaVOQBVfQLAE+2fnxGRBwG8EAAdACFDhnXXqotxrLrKZ1ioRRmoiKwBcC+Adap6PPhuM4DNADA5OfnLMzMzA58fIU2lqAEt2m6xV7LGKzL/rENaVR3iKoPaloGKyFkAbgfw7tD4A4Cq7gSwE2idAxjw9AghOYjp9NSBYTDOVVLpDkBEzgRwB4C9qnpz1vU8CEZIPYnlAXpdIYeOJCbTUHSsoo1nun1G3ajdDkBEBMAnADyYx/gTQvpLN+Gb0BAfPrpkE59rnDRjT/pHlSGgfw3grQC+LSL3t393tareWeGcCCEpZK1+165e1XMO4PDR4x1NX2JjJzmOrHnnKQttUpP4KquAvgZAqhqfENLiypv24uTpOcy3o8EHZ2Zx5U17cxnxNK3+EFv5m2E/ODOLjdv3YMXysSXGflnJliHPzqSJVJ4EJoTUn25X0N3Gy80ZmWxzljPKmhfVP+PQARDSUMxY2up75fgYTp6ew7rJicIGsmhnr0NHZrFi+diCfIMd6jLCHUFeQimKWB+Ak6fnsGXXPjoB0AEQQnJQ9GBVXr2emBMC8lf1xE4jp91jOYoyGMZqoBA6AEIaSiip0E8xtZOnO+Ua5rVl9GPG2Fbwh48ez7VS97kFy18ArfcJHQSA2ugV1QE6AELIAlkNU/LWxifp9Vjox/DhGr+CL3qYbFA6RMBoCc3RARAywuSRc/ZiamXLOdvKP8SHeMIV/MbtewC0dgl5Vure2Zw8PRcVh+t3r99hhQ6AEAIAS4ywL8+MlXn6n5NWwyuWj3U4Aavu8Y4mDAPNB+IEg+gxUIRREpqjAyBkBKmDnLONaU7ADHuYiLXPK8eTzwPkyQP4+a8cH0t0HMNssMtmWdUTIIQssmXXvkrE1MxQWhXOvLZ+t3J8DNNTEx3O4uDM7JLQzPTUxMJ1YVjnxKk5rFg+hmWChWvM6Nvz7LoQcxp1EpgzwtDSMMIdACEjSJ4wRVLC1zh89HhiDD/Er+jDE7/GiuVjHUY/vG+ZLA3/AIvVQn5+Se8X7mZY7ZMOdwCE1ABb+Xujlbbq7cdOYcemi7F76wasHB/DyvExrF29quNQmK14p6cmFnYK9t3a1as6wi3+Z3vW7q0bloRk1q5ehWXSuuauay7H9NTEQtjHy0GYE6CkQ7lwB0DICJO3mbu/1lb9eZqupz0rdjgrtnL3p4Dt88nTcwvlon4nEdb6p7WDtGu48k+GDoCQGpC3siTJ4ALIdRI2DTOYsTBMOL4ZYtsJJJE2HwvrxMpQLS+wcryVO1ixfCx3kjpW3grQEcSgAyCkYcROx4YrcU+ecwIx6Ya0cI3fIaQRykuHuxEa996gAyCkRhSVXwbQkXBNknJOk2sOV+JAckK225O6aViOIO2wluU8wlV9jFGq0+83TAITUnOKJnyTNHY8XhcnvN5CLlba6ZO+Nh9LyFrS2nPlTXsXGrqcODW38Dl8J18mak1gijqWoo1hSCfcARAyBITCaOGqNkycWpgmpsrpdXmMrNO6IfZ7X86Z97Ru2IAmRvh+RVb1XPnnp9Km8EVhU3jSJMJ6+lioxNi4fc8Sg+qbqfgTsqH+f6zaxsYywlW+xfyB5ORzUvVQWMFjOw77nCYFnSZD4efkoSOoYVN4Qkg2fmXuD0SFidy01TTQaVRDAwx0SjCnNWv35JVXDp/RTbOXPHr/fj4AMiuUCB0AIbUkybDnDbMUMX5pz4xVDIWJ4NipXmBppU4S6yY7V+6xZHDW3Dw2H5aBZkMHQMiQsEyWGkcfIvKVO7ZzMEllHx4J4/dG2tmDWFI5a5cQzm96amJB9M3PJewd4DE5ijwN65N6EJBk6AAIqRnhQSYz7PNabDXrdw+9SCrHdH0sDOX7++Yp0fSOyZ5tz4r18o0Z8Sx9IjaAzw8dACEl0GuYoWgzljQDCXQmVcPksRnutKRynjEAdISFbEfhu375ew8dmY2GtEKBuNjfhU9YZ512ZmlofiqtAhKRWwBcAeApVV2XdT2rgEhdKdMBhMlMC52Eydkw/GPG0eZh5ZY+xh47gRt7rj0nZrRj0gwW1vFGOtw5+BCVXb9764YluQTvmOy7MBTk/07CKiCu+JdS1yqgTwL4KIBdFc+DkK7otT9saMjtOV4JM0Ys/BI7ALZi+ViHIY1JNafp68Sqi+a1Zej9WKFgm30XOq/YqjzmfMK/v6S8BemNSh2Aqt4rImuqnAMhVROLac9rcjes0OmE5ZyhtIM//JVEeDjr4MwslsnSGn1/PqAXaebYO/mWlFklpezxWw5V7wAyEZHNADYDwOTkZMWzIU0lyciknVDNY5h8SMbCJ2nGOik+DuQvEY2R1o0rVOFMWsnbIbOV42NR6Wi7N4myW1KSbGrvAFR1J4CdQCsHUPF0CCmNcCVvK34z6El17L59o4V+rAeu7SbC5u1Z8gu20rdx/YnhcKUfKwtdJovzPXl6ruOwWuwdQvJKPcR+z5V/99TeARBSJXlj/LFGJEUPIlmyNhZaOXRktiOUE2ugDnT20LWdwpZd+zK1d6zz15U37cUyQUcSNzyNnLVKDw92Fa1wSoKhnvKhAyCkItL62hrmDGx1bvhqobWrV3UIuwGdydLDR493hJZiUs+Hjswu0ROy69PCUlaeGVYaGXlj+7F7Sf+p1AGIyGcBXALgHBF5DMC1qvqJKudEiCctNJHUnCTUounFqFm1TSjg5sfxTiLNyBtZukFAp1JoUmze7xT8PIxYxVE3O4FeK61IMlVXAf1eleMTMkiSVvoWetm4fQ/WTU501MPPa2eCNhbK8bX/G67bAwAd1Tx58PX54Wo/3F34e+z5edpDkvrBf2OE5CAtxj+INoVhpU9SC8e0lb9/lpeB9k4llGoA0kNAnljOI9wFdJMHYIev/kEHQEifiTVlARYNpo+9Hz56HMtkMSnrT8L6Z1hi14y4ddTK0zzdV+v4HUd4+Mw+7966obDAms1n7epVS/r6kvpAB0BIm7wrTFvpJzU8ybo/S8xsXhdF0sxw+lg7gAWp6FhVTt4DWjaO/+ydQJh4TsPu847LN5wvQ5eHzqN82BOYkD6zY9PF2LHpYqwcH8PK8TGsm5zAusnFfrt3XXP5wrW7t27A9NTEQqgkXDWvHB/LlIk4eXpu4RRvFqEzioWPtuzah91bNyzMzXoE21xsx7Jieesswsbte3BwZnbBOZmjpAGvH9wBkMbTbZVJ0mnX2POT9H7mtVMHP9baMXaI6+TpOUxPTURDPXadOQD7p63w/XOS8gW+2iiGF4tL+x2pN3QApFF0k0hMq94pIzGZFRKKVf6E4ZsYdo93COsmJxbuiyV8jROn5hacRax+3z6H1Uahg7FSUa7+6wlDQKTxWIjGwhv2uQzSNPVtRb5i+dhClU8s7h6erDVOnm5JNux9/2IIyUIysfCPOQ2bizWZScJ/x1X9aMIdAGkE3YR5sqp3vDRD2mneNGJ1/iGHjsxG6/NNljk2lq/u8Vh+IO1wV8wphLLSvqlMklIoV/71hw6AjDxJPW1D+mGs8vap9W0fQ8xJJBnnsJIHSD4A5u+3OH+enIC9Q3iqNyt8ReoNHQBpBD5BaavZJNJi/uHnpO/C9oY2h7Dbl5G0YvckOQegM9Gbh1if3zRCQ2/vmaYBROoPHQAZWWJaNEDL+NohpW6SwVnEdhu+I5Yd9gpj7HlO8SaR9z7fzjHJCcQSueHz0yQgylL/JP2HSWDSSExkbcuufR3tEq05uVW+eF17f7AJaJVs2u9C52B18sDiYS27pkzDaPH8LB0e393LnNHurRs6avmNok7IDqoBcfVP7gbqC3cAZCRJUuYEkg1wUp7Aeu/6kkg7ietDI/46oGVUY89M64GbpPmftDuYVywZM2y8bljzGBvH/h7CE8B5iJ1BKEv9kwwOOgDSKEJ5ZSO2OvchG298veHzipihkbbVNrB09xBiOw+Pn2NeQTYv8AZ06vh4Tp6e63AcoZS0vVtRfHWRVQWxEqi+0AGQkSSrHDOmn2NhC1s524q+iCaO4Ru6J82h16SpORxvvC22752P4U/3xiSlfbLaOz9v1KenJpa8V1KDdxr/+kMHQEaCvKdzw7JMr1RpRs4bzlBd04dVwlBNbJWepyNWbIUektb+0RNKToSlm3aNx4eqyjLaNP7DAR0AaQThAbAYtroN7zHjuEw6Nf8Pzsx2dMXyjsIkF2I5AEs+l2Ug00JDSXX6PlxlO5zQaGc1dWdbx+GHDoAMNbHTuv50bnhYySQNisS3fdhk4/Y9WLF8DLu3bkhczdv1YbglbOpimFOxbl4x8paI5lH3NPx8fOUTQCPeFOgASGWUIcyWha3OrflJeBDLh2fC+LaFdsKTsmZkY3OwBu1h4tjG9LmH8JBYFnmcQF7HZge48vYOoEMYTegAyFCS53SuGf6i/XH9vTFMwtknS8PvfVjGG27/syWYrTKpbPa+//IlpalAPD6f17FyhzBa0AGQgVOGMFvSwafwFGrWijipXaFfkcccx4lTcwvhJG+8s5K1sV1B0n1pz8mDl2swh1ZUBoOMNnQAZChJU+FMC6ck1bhbYtbIu1uIVdnkJUtuwTdv70Z0zRt9G+/EqbnUhG5Skrzbpjmk3tABkIGTJLKW556N21uJ0tBYhgYq7eSs59CR2US9/SzCUE9ekqSd/ffA4u4kaQw7aJV2ja/w8U4trESigW8miQ5ARO4E8E5VfaRfg4vIZQD+AsAZAP6rqt7Qr7HIaOANGpC+cs2bNPUxcksET08tTRgPClvt+/LT2Hv4XYHp/MSE7wx/MrfowbO8TpuOY7hI2wF8EsDdInIrgJtU9SdlDiwiZwD4KwC/BuAxAN8QkS+p6nfKHIfUlyI7gVDbJxRh85U7afH7rNV3nnaJvZLlUMyQp638DZtvzPjHwka28k9a6bOBe7NIdACqepuI7AHwAQD7ReRTAObd9zf3OPZFAB5W1e8BgIh8DsAbAdABkCixBuQxYqWVtorOG0u3nUY/Vv9eeC3U7skjO5Ek4RxiuYNw19QLWSt/hpCGi6z/2n4C4ASAcQBnwzmAEnghgEfd58cAvDK8SEQ2A9gMAJOTkyUOT6qmiNHYsenihYSrXwHHDoB5zZuQrBV9aFx70ehPwj8vz6o9D6Hj8DuYE6cWy0z9SWZg6S4sKblORpO0HMBlAG4G8CUAr1DVH5c8dkyAdsn/aqq6E8BOAFi/fn0fNuSk7tjpWjOcoWhbSCzuHWu44pOtSSeE+xECyiKcY+gkpqcmlszXi7nZ770zCQ172XST2CfVk7YD2AbgTar6QJ/GfgzA+e7ziwA83qexSA3pxWiYBo9/VugorDomdioXwMK14UGpcJxBOYFl0qrZt/mETdaBpUqd/ndA567KdP5j8w//rmnAm0laDuBX+jz2NwBcICIvBvB9AG8G8Pt9HpMMIRa2MAO/bnJiicHasmvfkgbnWeGUPBr7MePpQ0x5evnmZV7R4azssJnXEgr7CpgjsJLOpBPFg9Lmp+MYLiprCamqcwDeBWAvgAcB3NbH3QapMTs2Xdy14fD3hnHwPKeA81xnWKLVnIbtHnrBnmmtGdMqlGLs2HRxx07It2cEOpPZ4WG3pOfRiDeHSg+CqeqdAO6scg5keLCdQIgZNet+ZYli3yQltkovunL3xtiLqZ04Ndd1qMiv9m1nEeoGhaWtRloS3XZLSX0OCAF4EpiMIGF+oAixpGuMMG8QGv8kAbjw+aHjiP0cqoj6fr5J+INyjOuTJOgAyNCS1gsA6GxIbqEWH/ePGXufdE1b2Set9v1q3eZlY5pjCnvm+rkW2ZWkJW7LNPZ0IKMLHQAZGmKGKEtGOdTJyUoMh/IPWQ1VwtV8GqH4m6/Jj0k0eGfgw0O+pLOIjhIhIXQApHbkbSoeNi4H0HGflzwwfFgG6DTiyyTZQdh1fkcQXmv9g/3K/ODMbMd15mDmNXll7d/Hxrbnx8JT/TLwPN07+tABkJ4YhFEIDdHG7XsWkqOhNpBnx6aLF64NiSl5pq3203YEaSWWJtIW6xucRKxBTLhDoBEmZUAHQGqDrebN2B6cmV0QJ8uT1PUGPa2/rqdo9U7sejPYvhOZVfWYo7Fy0fAQV8yQm8OwZw2ifj8GD4eNPnQApDDWZtCHJGItGcsyGGYAbTXvq2MALEmuFiGP8feSE142okhf3SyRt5jzA5bmFWiESZnQAZBKSKtcMUM/PTXR8X2SsbY2h+aYej2ZG6u/jxn5cJzws4/923O86NrBmdnMg1nzio7rqtwJkNGDDoDkJlylWlLUG+Akrfm8mv9eCuHQkVYIyKplrEmLl07uR3gkdDTesHsRuiIyEjGJao/X4w8lLQjpF3QAZCCEDVxiTsIbd2BR0TK8x7BVtD8YVYY2j6ltdtMLwPccjp0zyIqrh8bfnA1X4aQfiOrwLDXWr1+v+/fvr3oajSerTDOUJADQcfjKG7kwzBIzuuE9wNLTtWldwHrBnIFPzKaNkyQUZ60mPeHfYVjS6qEDIL0gIgdUdX34+8rE4EgzOHl6Llpbb+JnwKJxy6r0CTtpheJs9qzpqQmsHB/rSN5mHdJKwmSVTUzt8NHjC/NdJq33sPGmpyawe+sGrJucWFAstXmEImuhiJsPrx2cme0Yi8af9AuGgEhhsgySzwdk1eCfODW3ZMcQCqKZobTqG1vxhzsGyxmYIqbPPRRR7UzS8UmL3yclc8Pkb0y2IinuT/E20m+4AyADwVbgMaNmq137ft3kxELIxIyircSBReexe+uGjpW+5QzCVfPho8dz9do1/LXrJicWVvo+Uet3BhbmAhZP/xat3PFjDkq7nxDmAEhfCFe8ZqRNshlYNORJ8XFf5x9T27Sm57HvAOCuay4HUKyyJhzLf7YqJD93e7dYDsLeyxvytMqovBIYhBQlKQfAEBAZCGHTEgAdVT/eIIY193mlHMLr8/S/nZ6aSDx0FY7tdymhomeSkS8KjT8ZJHQApGfSDnWlGUOL4x86Mtuh2dNtwjYkrL03DX+/qvfOxodhQueS1JQlbdwYaZLNNPxk0NABkL6RVtZov++m1j4voayCEQsZmayD4Us/Y2GZtHBN0eQttXZIVdABkK7JIxcc1rn7TlpJ9fq2Su8H/qCWzc9W7BaSWjc5ERVxyyLMbdCwk7pDB0BKJ8kxAPmVN8uSQvBnAYBOATmgM/cQk5aOVQ8V1fNJgnr7pGroAEjX+LJI/zlmEC2k0gthhU7e/r0xrILISzPHxsuzmg8rgyjfQIYFOgBSiFgZY5KRDI1rVjvGGGl6/UWMv2+/aIRJYluJZ0lLh44vvD8v1NsnVUMHQAbCusmJqNaPGfgwNh+jrLCQX7EfnJldqNf3Oj5+vKx8BA05GVYqcQAi8iYAHwTwSwAuUlWe7qo5afHqmOEL4+ShcQUWD3MByavoPE3ZV46PYffWDR2tIo1wdxL+HlisSPIJ4bQ5JO18uoFOg1RJVVIQhwBcCeDeisZvDKbJUzVm/E26YeX42EKZpTkJo2gFkMkxeKmIvOGh8CCXOaJwDl7eISlnQOE2MmxUsgNQ1QcBQKRPtX4kF0VWn2lhjtj9tjOwVbo3zrZCj9XL5zH+tnOwuP7G7Xui+YXdWzcs2ankad8IpO88ijRoT7qGFUCkDtQ+ByAimwFsBoDJycmKZzM81MXAJK3E/arbtHqAuApniNXo2zutm5yIqn1u3L5nSYP1mCaQNXT3Gj5h1y9W9pBRpG8OQES+AmB15KttqvrFvM9R1Z0AdgItMbiSptdo8jiHWEVP2jVJO4MkeefQkM7r0vp7O5TlD495/Encjdv3LPl+XjvlmNPkHcKdgZePiDWDCd81JOvvmIljUgf65gBU9dJ+PZsUJ2/ooyzCRKkRE2pLUta0Od91zeULp4iT5BnCKiMf2w9zIJaQTurGVZawGyF1p1I5aBG5B2BMNKIAAAzYSURBVMB781YBUQ46P2lyzOE1MYMXKmSGksihIqaXdo6pZnr8NUCnA7D4flbYxj/Lv4NVAk1PTSTG3Q3vHNLUPXuBDoTUgVrJQYvIbwH4SwDnAtgjIver6oaM20gBwsNKg9an8avvUBAuXJX7lbuFfrwDiom3JTkEE3VLS9KmvXs3Xbho5MmwUlUV0BcAfKGKsckiadU7thJPS4KGOYBDR2Y7+tr6+L+dBbCSzbAZTFJzeN/43XICwGJsvohxj+0IktRKy4JOgdSZ2lcBkd4IJRB6MUimnVPkeu9E1q5eFa3WCROuaZ3ALHQUVgKlhbH6tUKvS6UVId1CB0AWCA1amGy98qa9HYnYsJLFjLLPCfg8wOGjxzuuATqNpe+ta5hWvw8nxU7x9vKeAPvwkmZCB9AQujFstlq3OnkzvBbisdLLsJmKr+mPrfhDwoohf1DMVvw+XJOVsB1UiSVLOcmwQwdAFoipePq2jSEm4hY2fPFhHMP0evIYS3MgfoXuq4Z8JVI3RpeGm5AWdAAkmsz1XbvS5BnmtVMW2ap4/PfAYrzeh2ySDpsZPnQUNmvJY/wHZdjpQMiwQgdAlhCu4P3hq7SafF96mkcCwojV59vvwwNsPn/g7+1lJ5AGdwlklKEDGFGKCJWFVSwWqvFhFm+kfZLXSzsknf61e4wrb9q7JJZv44bESkNZq09IOdABkAWyKmr8d2a88zZpOXl6LrEtZNIhtbBCKKwI8gJuZRt4lniSJkAHUHOKGp60huzhc9J0cPx3oVKm/TN2j5eA9lhuoNsevknNWrLo1pAPWjuJkCqgAxhxslb1Pulrp3ezEqxhS8U0o+oPc/ndQqwBTOz+WM3+ll37+l7JE/6dceVPRhE6gJrS7co1KYGa9pxY2Wb4PKv0iQmoJY1vhPX7QLqej1GGVENRR5EkokfIKML/ukeAWDllrHTSfg7LML3x9/H40FiaFEQRoxpKNIQni4F48tfI6r7Vzcq8yK6hnzpBhFQNHUBNKbpytdO6dl1sxR1KOOQhdBLhOFnzD8kzfmz3U9ZOIO91TPqSJkAHMMSEhjI0zvZ9Vsw+LVHs77dnHDqymCsoQvjctJV/SNr5gV4T5DTypKnQAdScItUqJ07NdTgBb/h7GT8MKYWlnFmtIYtQ5rN6gU6BNAE6gCEmZpz9qtx+DhOasTi63wWEIR4vEeETxr5pS6y+354VG7MXykiQlz0nQoYROoAhJ+y6lVTnH1LUCPqkbYhJPHsJiCJ9A/z1DM8QMjjoACqiTANnxjnr8FLeOvuk+fnmMjFdIC8CF4afyn7fXp5Jp0JICzqAESEtIZvWLatInXu4qg9bRsbmEDsrkDZHrvwJGRx0AAOm7FBHnueljWESC16+ISlXYKTV9idVGPXDoNNJENIbdAAkkywnE9MQIoTUHzqAAVP2yjjteXkPVGWdts0zfjhmLw4h72EzQkhv0AHUkF6cQ3hvmBjuxth347R6PaRFJ0BI/6nEAYjIDgC/AeA0gMMA/p2qHqtiLlVR1DAWrXG33yWViPYyp7R5+rMBRUk71EYIKR9RzdnRo8xBRV4P4KuqOiciNwKAqr4v677169fr/v37+z6/GIOoTgmNaJHm57HqHl+maVU6/Vi9x8Y3rH9A3nH9obYi9xJCkhGRA6q6Pvx9JTsAVb3bfbwPwO9UMY8sqi5JNGOYpbmfh36rWtq8Nm7fA2BR6rlIY5WiOxZCSG/UIQfwdgCfT/pSRDYD2AwAk5OTg5rTAhaGGMQJ1VDOOVT07PberLJOo4wSVZOE6LZ/byhLQQjpH31zACLyFQCrI19tU9Uvtq/ZBmAOwKeTnqOqOwHsBFohoD5MdQmxhGQvz+nGWWQlawe5Oykylp0W9g1kisKVPyGDoW8OQFUvTfteRN4G4AoAr9MqEhEFsMNSgwpL9DJG3gNgafdmyUXngY1UCKk/VSWBLwNwM4DXqurTee8bdBK4W2niWDI3771lPrvbXUJo8PvxHoSQwVGrJDCAjwIYB/BlEQGA+1T1qormkouqjF4voZ5u59zLyp8QMjxUsgPolirLQLuh28NQecI4dc0BEELqR912ACQD6uMTQvoNdwA1IO0AGOPvhJBe4Q5gyOhFoI0QQvJAB1AD2AyFEFIFdAA1h86AENIv6ABqQJUrf+46CGkuy6qeAOmdLbv2sWafEFIY7gAqpMpST5aZEkLoAIaYLCNOo04ISYMOoELKNtSxfr+DGpsQMnzQAQwxsR4A/uwAwzuEkDToAAZMzBj32os37BzWzU6AENI86ABKpFv56F4JO4fxFDEhJA90AAOiX2EZxvIJId1CB1ACoXG/8qa9Cz1x0+rzrd9wv4w2nQEhJA06gAERW6mXeXiLxp4QUhQ6gBJIM+4xw2wrf1bpEEKqhFIQA2bHpovZMJ0QUgvYEKZCuPInhAyCpIYw3AEQQkhDYQ6gQrjyJ4RUCXcAIwjloQkheaADIISQhlJJCEhErgPwRgDzAJ4C8Aeq+ngVcxklKAJHCClCVTuAHao6raovA3AHgA/0e0CGRQghpJNKdgCqetx9XAlgeGpRawx1gQghRaisCkhErgewCcA/A/g3KddtBrAZACYnJwuPw7AIIYTE6ZsDEJGvAFgd+Wqbqn5RVbcB2CYifwbgXQCujT1HVXcC2Am0DoL1a76jBJ0bISQPlZ8EFpEpAHtUdV3Wtb2cBObKnxDSVGp1ElhELnAf3wDgoSrmQQghTaaqHMANInIhWmWgMwCu6veAXPkTQkgnVVUB/XYV4xJCCFmEJ4EJIaSh0AEQQkhDoQMghJCGQgdACCENhQ6AEEIaCh0AIYQ0lMpPAhdBRJ5G69xA2ZwD4Ad9eO6gGYX3GIV3AEbjPUbhHYDReI9e32FKVc8NfzlUDqBfiMj+2DHpYWMU3mMU3gEYjfcYhXcARuM9+vUODAERQkhDoQMghJCGQgfQYmfVEyiJUXiPUXgHYDTeYxTeARiN9+jLOzAHQAghDYU7AEIIaSh0AIQQ0lDoANqIyHUiclBE7heRu0XkvKrnVBQR2SEiD7Xf4wsi8ryq59QNIvImEXlAROZFZKjK90TkMhH5rog8LCL/ser5dIOI3CIiT4nIoarn0i0icr6I/E8RebD939KfVD2nbhCR54jIP4rIt9rv8aFSn88cQAsRWaWqx9s//wcAL1XVvjeqKRMReT2Ar6rqnIjcCACq+r6Kp1UYEfkltJoFfRzAe1W1uz6gA0ZEzgDwfwD8GoDHAHwDwO+p6ncqnVhBROQ1AJ4FsCtPq9Y6IiIvAPACVf2miJwN4ACA3xzCfxcCYKWqPisiZwL4GoA/UdX7yng+dwBtzPi3WQlg6Dyjqt6tqnPtj/cBeFGV8+kWVX1QVb9b9Ty64CIAD6vq91T1NIDPAXhjxXMqjKreC2C26nn0gqo+oarfbP/8DIAHAbyw2lkVR1s82/54ZvtPabaJDsAhIteLyKMA3gLgA1XPp0feDuCuqifRMF4I4FH3+TEModEZNURkDYCXA/h6tTPpDhE5Q0TuB/AUgC+ramnv0SgHICJfEZFDkT9vBABV3aaq5wP4NIB3VTvbOFnv0L5mG4A5tN6jluR5jyFEIr8bup3kKCEiZwG4HcC7g13+0KCqP1XVl6G1o79IREoLy1XVFL4SVPXSnJd+BsAeANf2cTpdkfUOIvI2AFcAeJ3WOMFT4N/FMPEYgPPd5xcBeLyiuTSedsz8dgCfVtXdVc+nV1T1mIjcA+AyAKUk6Bu1A0hDRC5wH98A4KGq5tItInIZgPcBeIOq/rjq+TSQbwC4QEReLCLLAbwZwJcqnlMjaSdPPwHgQVW9uer5dIuInGvVfCKyAsClKNE2sQqojYjcDuBCtKpPZgBcparfr3ZWxRCRhwGMA/hh+1f3DVslEwCIyG8B+EsA5wI4BuB+Vd1Q7azyISK/DuAjAM4AcIuqXl/xlAojIp8FcAlaEsRPArhWVT9R6aQKIiKvBvC/AHwbrf+nAeBqVb2zulkVR0SmAdyK1n9PywDcpqofLu35dACEENJMGAIihJCGQgdACCENhQ6AEEIaCh0AIYQ0FDoAQghpKHQAhHRJW3Hyn0Rkov35+e3PU1XPjZA80AEQ0iWq+iiAjwG4of2rGwDsVNWZ6mZFSH54DoCQHmjLDRwAcAuAdwB4eVsJlJDa0ygtIELKRlV/IiJbAPwdgNfT+JNhgiEgQnpnI4AnAAxl8xTSXOgACOkBEXkZWh3AXgXgT9udqAgZCugACOmStuLkx9DSmj8CYAeA/1TtrAjJDx0AId3zDgBHVPXL7c//GcAvishrK5wTIblhFRAhhDQU7gAIIaSh0AEQQkhDoQMghJCGQgdACCENhQ6AEEIaCh0AIYQ0FDoAQghpKP8fIPVHspjC7+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.mine import MINE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100       # batch size of data sample\n",
    "lr = 5e-5              # learning rate\n",
    "ma_rate = 0.1          # rate of moving average in the gradient estimate \n",
    "\n",
    "mine_list = []\n",
    "for i in range(rep):\n",
    "    mine_list.append(MINE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                          batch_size=batch_size,lr=lr,ma_rate=ma_rate))\n",
    "dXY_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    mine_state_list = checkpoint['mine_state_list']\n",
    "    for i in range(rep):\n",
    "        mine_list[i].load_state_dict(mine_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEICAYAAACtXxSQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5gT1frA8e+7C0tHqvQqooIi4AIqoqKI/WK5KnaxYEHkggpIlSqgIkUs6M+GXa9eFRsKqFhAuiBIld5Bet3d8/tjZpfsbpJNmckk2ffzPDxMJjPnvJls3pycOXNGjDEopZRKbCleB6CUUip6msyVUioJaDJXSqkkoMlcKaWSgCZzpZRKAprMlVIqCWgyTyIi8oOI3OtxDKeIyHwR2Scij4Sw/ZMi8ra9XFtE9otIqvuRJi8RuVVEpngdh4otTeYJRkTWiMghO+ltFZHXRaR0mGXUFREjIkVcCLEn8IMxpowxZlw4Oxpj1hljShtjMl2IKyn5ey+NMe8YY9q7VJ/nDQblnybzxHS1MaY00BxoAfTzOB5fdYA/vQ7Cl0tfWkrFFU3mCcwYsxH4Gjg973MikiIi/URkrYhsE5G3ROQE++mf7P932y38c0SkgYj8KCJ7RGSHiHwQqF4R+ZeI/Ckiu+2W2mn2+mlAW+B5u9yGfvatZ9ezT0S+Ayr5PJfTyhSRjiIyJ8++3UXkc3u5mIg8IyLr7F8oL4lICfu5C0Vkg4j0EpEtwOv2+p4isllENonIvXZdDcIo71H7WG4WkU4+cZUQkWftY71HRH722fdsEfnVPlYLReTCIMe1uoj8V0S2i8jfvt1UItJSROaIyF47vtFB3su7RORnn32NiDwkIivs4z5ERE4Skd/s8j4UkTR72/IiMtmO4R97uab93DCgjc/7+7y9/lQR+U5EdonIMhG50afuK0RkiV3vRhF5LNDrV1Eyxui/BPoHrAHa2cu1sFrBQ+zHPwD32st3AyuB+kBp4BNgkv1cXcAARXzKfQ/oi/UFXxw4L0D9DYEDwCVAUaxulZVAWt4YAuz/GzAaKAacD+wD3s4bF1DSfu5kn31nAx3t5THA50AFoAzwBfCU/dyFQAYw0q6nBHAZsAVobJc9ya6rQRjlDbZf8xXAQaC8/fwE+3XXAFKBc+16awA77e1T7GO2E6js57ikAHOBAUCa/b6tBi71OW6328ulgbODvJd3AT/7PDb2aytrv/4jwFS7jhOAJcCd9rYVgevtY1QG+Aj4n09Zud5foBSwHuhkv2/NgR1AY/v5zUAbe7k80Nzrz1Cy/vM8AP0X5htmJfP9wG5gLfACUMJ+LueDZn9YH/LZ7xTgmP2B85cA3gImAjULqL8/8KHP4xRgI3Bh3hj87FvbToqlfNa9i59kbj9+GxhgL5+MldxLAoL1hXKSTznnAH/byxcCR4HiPs+/hp2c7ccN7LoahFjeoTzHaxtwtv36DwFn+nm9vbC/QH3WfZudOPOsbwWsy7PuCeB1e/knYBBQKc82/t7Lu8ifzFv7PJ4L9PJ5/CwwJsB71hT4x+dxrvcXuAmYkWefl4GB9vI64H6grNefnWT/p90siekaY0w5Y0wdY8xDxphDfrapjpXss63FSuRVApTZEyup/W53odwdYLtc5RpjsrBaZjVCiLs6VmI4kCeuQN4FbraXb8FqIR4EKmMl9bl298Vu4Bt7fbbtxpjDeepe7/PYdzmU8nYaYzJ8Hh/EaiFXwvols8pP/HWAG7LLtMs9D6gWYNvqebbtw/H36x6sX0V/ichsEbnKTxnBbPVZPuTncWkAESkpIi/bXUZ7sb5EykngEUZ1gFZ54r4VqGo/fz3WL5O1dvfaOWHGrUKkJ4aS1yasD1q27FbxVvwkXmPMFuA+ABE5D/heRH4yxqz0U+4Z2Q9ERLC6ezaGENNmoLyIlPJJ6LWxWo7+TAEqiUhTrKTe3V6/AysBNTbWeQN/8pa5Gajp87iWz3Io5QWyAzgMnAQszPPceqyW+X0hlLMe65fAyf6eNMasAG4WkRTgOuBjEalI4GMXqUexfsW1MsZssY/9fKwvevzUtx740RhzSYC4ZwMdRKQo8DDwIbmPvXKItsyT13tAd7FOOJYGhgMf2K3L7UAWVp8pACJyQ/aJLuAfrA+tvyGCHwJXisjF9gf0Uaw+2F8LCsgYsxaYAwwSkTT7S+PqINtnAB8DT2P1ZX9nr88CXgGeE5ET7fhriMilQar/EOgkIqeJSEmsvunseiIpz3ff14DR9gnMVPskZDGsbqKrReRSe31x+2RqTT9F/Q7sFeukbQl7+9NFpIUdz20iUtmub7e9TyZ+3ssolcH6YtstIhWAgXme35qnrslAQxG5XUSK2v9a2Mc5Tawx7ycYY44Be/H/N6UcoMk8eb2GdZLvJ+BvrNZjVwC7q2IY8Iv90/hsrCGOs0RkP9bJsm7GmL/zFmqMWQbcBozHapVejTVU8miIcd2C1T+8CytRvFXA9u8C7YCP8nRz9MI68TrT7g74HqtF6Zcx5mtgHDDd3u83+6kjkZSXx2PAIqwTtLuwTrymGGPWAx2wuku2Y7ViH8fP585YY+uvxuqj/hvr2L6KdYISrBO4f9rvz1isE8GHA7yX0RiDdcJ4BzATq7vJ11jg3/ZIl3HGmH1Ae6Aj1q+2LRw/8QxwO7DGPqYPYP3tKBeIMXpzClX4iDWccjFQLM+XhFIJSVvmqtAQkWvtn/7lsVqPX2giV8lCk7kqTO7H6u5YhdV3+6C34SjlHO1mUUqpJKAtc6WUSgKejDOvVKmSqVu3rhdVK6VUwpo7d+4OY0xlf895kszr1q3LnDlzCt5QKaVUDhEJeMW0drMopVQS0GSulFJJQJO5UkolAU3mSimVBDSZK6VUEtBkrpRSSUCTuVJKJQFN5kqphDBgwiDGPJN3evX4NHXDHDKyYjuHmyZzpVRCmNioAyPOutbrMHK0/P4zmkz9Nt/6qRvmcOuKIgxcODmm8eht45RSKgzPf/46P6Rksq5Uut/nNxzcA1RkxaGsmMalyVwppcIwtEyzoM+niHW71OwJaTOyDCKQKhJkr+hpN4tSilkzpjJ01BOulT9uzCDefG2MK2W3++Zdbv30BQC2bN1Ehy9f4+kJQ8Iq44X/e5qq0xcw8tURUccj9r2vjf1/zR8X0mbWX1GXWxBN5kopJmxbyvMtbmLYyF6ulD/8zA4MqtvClbIXF2vE1HLnAvDR/95iVsnmTDq1ZVhlLCppdVIsrHRCAVsWLLsB7nuniNWHjvjd1kmazJVSrClzIgBHShZ3rY6DUsq1sgFGjR/CsIaXOVrmT1OnMOjN59i4cROzZs8MaZ+clnmM7/ujfeZKqaQwo37NKPb2n3m7ksXW2m2ZvHQO61Nrs6WAUvYfPsD2QweAcgFKdI8jLXMRWSMii0RkgYjoROVKqYAeevdZhj3d19U6JE8qPXjgAD1fG8GixQsCbG/Jm4C3plQFYH1qbQD+8+H4oPU2+G0FT+2o4bcstznZzdLWGNPUGON/vI5SSgGfVLuY8ek3cPdH41ixfCn//vwVur01yvF6+r80hMETBgPw7KQxvFXvMgavngXA5o3rueqrN6gzLXDXyeA3n8u37v3KbUKuf+axeuw+vDfMqCOn3SxKqRzbHTgBGKqvKp1PpR8/4ecGV0IZGBvm/j3eGMmsGg2gyEn5ntuaUpVXTrkagAHA4VSr7b2kVH2qTl9AtaxNbC7RNN9+BmHLpo1UrV6DVWXSQorjwMEDzF2xkPPPPDffc6f+tjqMVxQdp1rmBpgiInNFpLO/DUSks4jMEZE527dvd6hapZSTPq16EQAL582m34RBbNqwIaz9v/rsPapOX8Bjrz0V0vZZqflT0JDnBnDKtB94cdxT9J8wKOC+79a5lFV+Erk/R1NTAdiZYt0+c3NK9ZznLvt6EgfTrMT9wwln03TZdlYuXxJSuQD3/vgBN+4qSd1pv4W8jxucSuatjTHNgcuBLiJyft4NjDETjTHpxpj0ypX93o9UKeWRvJezvLj8Z15t1IHnv3wjrHIWL/sDgG/qnAXA8Kf78Pj/hZbYs/3Q6DT2SDmGn34xrzTqwKhn+4e1vz/zq9cK+NyC4mfwTYXzcq17Yfa3fFu+dUhlTy/eHIDDUiLyAB3gSDeLMWaT/f82EfkUaAn85ETZSqnY21nSGka4r2SxkPf54O2XOSZW+3BHijXUcVz6jbm2mfr9R5B6st/9xz03kK9PPYm/ijUC4JhYreXDqcKI5wZA0+us7cYP42DmMTjzmlz7rytezW+5b73zAlIxvO6jd6u3DWv7ULT6bQmda1XmnpruNGajTuYiUgpIMcbss5fbA4OjjkypJPL15++zaOFcevZ/2utQ/Mo78iL78Rc1WlNnZB8yyKR3r5G5thk5ojf7yxRnSJcnAehWoxXUaBW0nlsDJHKAb06tz/ziTfKtP1wijRcaXpXzePjpV/rdP3vkSV6TKp7AomKNg8YVC2sPH6Xvio10qlEp55J/JznRzVIF+FlEFgK/A18aY75xoFylksYrsp/R593Kj1O/ik19zz9F1ekL6PXKMACGPNOHJlO/ZfyzuaeQHTryCaZNCRzTYSnJMy1vZEzLm/M9N6blDbzS6Bo/e4Vuct3jyX9e8TP9bvOaTyKPRDwkcl+L9x9ypdyoW+bGmNWA/3dBKQXArFLWyImdO7a5VsfIYY/z3Lm3cteKyZQ+eACAqfWslu6M0xuxLaUKy6uUYfnSBTQ8rSkjhvfk+XNuYdHeWVCifIHlz505g7TiJTijqTX62EhqznPnT/kQijbMtX3V6f7HdPvaLRVCfn3JYuuRY1DG+XL1cn6lEtSgMf0ZO2pAzuP95ax+7pl1jndlZPokXICParTjss1Woj9WrCgAs8o0YblPIh70nP8Tjvcc2M8l//hv/y3Pk8hVYAez3JkaV5O5Uglg/LMDGD70sZzHLz47hBfPvJ5Pm/p0Idgd3VkivJxu3cTBdwhetuw5UrL7xfOOwnix6fX8XCb/RFVbUo6fYHx94rM5ywPHDsi3rYo9vWhIqTg1YPxAJp/Wku6r/2LSmWexLrUOp3/wBktWLGbCudYokfVFauRsn31KLfRWcuQn4b6odnx0yOTG4c1QqNyhyVypGDAFJM6rvnqDTWmVGZeSxnkXXcJ/3hzF+6dbrevHTz6esDuf2BROzH/l4tJFc5Ewp+nbUq3gfnJ/nhzbn82Nmuc8zhL9gR8Ot2ZT1HdBqSg0mDaDFt9/HnSbwc/2Icvuuy5Rwv8Us3NKNGVTag3uJZURgx/j/drtQ6r/gJTmiZeG0HZHKm82vjTgdhuLVsn1ePzT/fik2sUh1ZHXS02u5+8i9XIe++vKUbGnLXOlgGHDe/Juq/bc+suXjGtzJydlrOKXS64vcL/9Uob9qWUYMag7m+pXp8qarfTt/0zO808N7cELre8IOY7dUoHJZ58dVuxvN7SS+GEp6ff5K796k50lcg84G5b+77DqUPFPk7lSwJLT6rIzpRLrGlpTnRY058eoJx9l1ak1oYp1peCGBjX4uHo7qAm+k7tmFg08WdOg5/qzrmYl/u+GbrnWryzaIKzYs6+UDGRuCR05XBhoMlfKh7/uzJefG0xakTQ6de2ds278+TflSqLHUo8PAaw6fQHnHJjLp1fdE/Qc44tNC275q+TjzsBE7TNXCa7H6yMYPqhH0G2GjujJ3R+P44/fZ0dUx8Cm/2JE49xdH3lbw2tPqJTr8W+lrImmyPKfzZ8e1jNnecyIPhHFpRLTr//sd6VcbZmrhDViUDfePb8TdWuuIVg6nNTyCvZIORp8/RZNWga/qfC20scvzXt2SE+aND8HStZjj5QDYNSgHrzSpgNI7ombFhQ/w295+yqWzvV4/h9z+PjIVr4895acdVPOPCVoTCq5LDtw2JVyNZmrhJV9OfmelNzXRnd7axQnbvmHPxvWpvGSv9lzjpU4x51/h9+k//A7zzC1ejsAfit1/EZZT593C232/Q5YIzc+mvQyo88P/WTm15+9z5sNck8KNa71bfm2CzQniVLh0GSuElag4bof1GoP9vTVy1v6v8lvjzdGsr5CeZouWcXHrfJPIpVtVunjibZrzeAzAubVqeypYW2vVDQ0mauklvdinVrTfueupV/xrj3b34xWevWiSg56AlQlrOw0/U9KRYaOeMLvNhmSu71yTNLCmrb1qIR+cwalvKTJXCUs326WmU38n0QMdMMCpbxiAnYQRkeTuYpLn384KaztD6dYQwUHPP+kC9Eo5RwdZ64KjaGjetO58hk88dKQkPdZXKwRDabNYGLj6O58o5Tb3JpoS0+AKs/0mzCIvWVLMu72x3Otn3qmNSvgshrVGPJMH8ocBkMGG2pU4LaGF9Cs9Tl+y9svLty+RSmHtbJvIuI0TebKVUOHPMamk6pxyvpd7D+wn6OVy5G2fQ99B4/h1UYdABiXZ5+ladaQvo3FK/PrWem5nkv562umTfsvh9KK8Hyb0Md8KxUv6pVw56S6JnPlqvfPbc+OlBOpXGUb21NOBKBW5vpck1ENH9yD0sVK8EivYbn2XVukTr7yJtW/nNp11lHvwGY3w1bKNS71smifuXLXQXta1uxEDrA+tRbDBzyS83hcmzsY3vKGkG4ADLAutTY/lUkveEOlChFN5iosVacv4LZPXihwu+FDH7MnwPI/0dS4tnf7XT90VG+/6/MyeW5UrFSicKtlrt0sKmzflz+3wG38zUESiudbdIxoP6UShVvNEMda5iKSKiLzRWSyU2WqxNRvwiCvQ1AqblVMc6cN7WSp3YClQFkHy1QJ4vfpPzJ29xIaLV/Dqy0DT1ylVGFXLa2oK+U60jIXkZrAlcCrTpSnYmfEoG70emWY3+eG93+Ydt+8S5+Xh1qPfU5aZrvs60nUnfYbHU0qU8udw3vp7VyNVynln1Mt8zFATyDgVRsi0hnoDFC7dm2HqlWRGNavC4eqVWbwff35sHV7NqXWoPTTffi6aTpX/jiNvkOfB+DXc9JZXKwRixs2ovTgHky6oEO+svLelGFHSuWYvAalEtWcvQdofoLzFw5FncxF5CpgmzFmrohcGGg7Y8xEYCJAenq6Wyd0VQhmtD6bBcXPoMSonuw89wYAfmhyBquL1GdPnb8A+O37acwp0TRnn3F5LtCpOn0B9TNWQ5H6sQtcqSSwL8Od2Vmc6GZpDfxLRNYA7wMXicjbDpSrXLI5zW49Fy2C2AOlxJ4wwqRYQwnnL5hRYDmrNZErFTYJcpPvaESdzI0xTxhjahpj6gIdgWnGmMjGpamY8L1hwxGs2Qazk7px6Q9NKWVx6yOmFw0VQln22/5hetuci2/2p1pXamaJ8Ot3Uxl8Vv7+caVU9FJcSueOJnNjzA/GmKucLFM5LzuZ+9644e8i1k2LM1NT+P6PqZ7EpVRhkBKv3Swq8WRJ4Lc9KyWFmaf5v2uPUip6bnWz6OX8hcTP33zHm3v+JDM1hd0Vzw+43Uc1dJy4Um5KcekMqCbzQuKHBd/xRatbvA5DqULPre4Q7WZJMk8NfZxPP3oz5/FnH0+i/Tdvs716JQ+jUkpla1iquCvlass8iQx/qifjWt/KFTt+YuXAR3n2wtuhonWF5h+1Tvc4OqUUwEUV3Zm+SpN5kqg37VfKtbwQgJUnVGV16yreBqSUiilN5gnsvg/GsKlMeb684k4OSUkO2Xf1WV60oceRKaViTZN5AvvixAu9DkEpFSf0BGgSaDL1W69DUEp5TFvmCeis7ydT6/BWKHUWANtStH9cqcJOk3kC2phak42lanodhlIqjmg3SwK445MJPPTus16HoZSKY9oyj3NvvvwcUxq2BSDr/efYXaIElD3b46iUUvFGk3mcG93gtJzl/1Vp62EkSql4pt0scc53mlqllApEk7lSSiUBTeZKKZUENJkrpVQS0BOgceLDd/+PxTvXMbjrIACG9unCuiYNQE96KqVCoC3zODFd9jDx9Gvp++JgAJ6/5D4+10SulAqRJvM4sb1UaQB2lC9D74lDPY5GKZVotJslTvxcpiUAn1VpCzrVigpDw2PLddpjpcncSx+8/SpzD23lrQZXeh2KSmD1921jeQVN5oVd1MlcRIoDPwHF7PI+NsYMjLbcwmBqkX18rolcRSklK8vrEFQccKJlfgS4yBizX0SKAj+LyNfGmJkOlJ3UNpcp53UISqkkEfUJUGPZbz8sav8z0ZarVDy5f9GnXocQmHgdgIoHjoxmEZFUEVkAbAO+M8bM8rNNZxGZIyJztm/f7kS1CW9fkZJeh6Bc1mZfvo+CipEGx1a6Wv7le39ztfxwOZLMjTGZxpimQE2gpYic7mebicaYdGNMeuXKlZ2oNqENebYvf6Wd4nUYKoBSZj8nZaz2OgwVobtWfcW4yjVcreP1Dg+6Wn64HB1nbozZDfwAXOZkucloQvMbvA5BBVHKHOCXS66jbsYaa4XRnsNQlDF7vQ4BgJJHMmje7JyI959cci8VsxKrByHqZC4ilUWknL1cAmgH/BVtucms+xujvA4hpi7eHV8/R8MxrkRJJh38m+L7DnlS/91/feFJvZE6+cjfXodgifI8QulSZUklsUYJOTGapRrwpoikYn05fGiMmexAuUnpkUlP82Gd9l6HoULU8jzrvbrkymsZO31BBCVEl1VK7d1f8EYJ7vpNU/lv9Yu9DiPhOTGa5Q9jTDNjTBNjzOnGmMFOBJaMho7szYc1L/E6DA9oF0U8K2YOx7zOsw/MzVmecOujrtXTd/m3rpUdb3Rulhj6unlLr0PwRIljx7wOISaaHf7Dk3pNlK3/h2Y6M+yyXoZzXSw55yqi1PX+Xo6Ukwg0mcdI8++/ZFWR+l6H4Yn6f64K+Nw1W6fHMBJ3Vdu/my1tm0a8vxhv+mh79RkZ8zoL+vppukNHEoVLk3mMbEp1d5hUPDNGr2oJRQOHhkJWztrmSDm+TszaGvG+V2/7ocBt7lr5Zb51126ZFnGdofbsXbz71yBFJNbfrSZzlVC8ar36SjGZ0SWaAE7ZvZFb1n4T9n6SJ3M1OLQOgCpZWxyJK1on7ip4uGKPK+/I9ViA4kcyXIrouMv2ZfpdnxnFfDcdWRjxvtHQZB4DD7432usQlIMEw4s39whp25O2hDdWuYwjQyCtFmVKXA+tK7jpnKg/6Ma0vZMixP48kSbzGPi06kVehxCUGP+tk2RSM3O942XWCqHMokePcd+S/0VchxwLPyHHauyQkyc849EJmfFxAVSoNJm7bMgzfbwOwTUNj63wOoSQlcvY53iZlY/tCmm72y+8gZvWT6HFwflh19Gn39Nh72MkNk3apjtyJ/NUkxlyP7MU8I0jxhS4TTCh7mp8ruzN23U2plQVOm1InBP0msxd9P2nXzDhrBtjVl+HrdO5b8n/uGXNN3SZ91HI+0X60U/LCu2nZN4+XV/hfF7POTCXE9gdxh75nXbUwYuTQwy+YePGjL2jJ0WzCv4FFEn+8upEnfgkwgv3zuSaefP8bhdKfEWL5Z50zl8ir5a1KbwAw1DW7ObFm3tQ3Bzv5kpv3Yanbu/uWp1O02TuotvK1YppfcWOZjCky5OM7tSb/o8Oy/XcIz+9FXC/YMk2HG7PUgdQNSO6+TKceq3B3LIm/JOYTsp+hXUObXa3Hp9fAO93eIDHeg93rOwbqzTIt67j4t9D3j9Bu9ujosm8kOgz0PmTsL6J8YE/PuHyBXODbB2eR6a/5lhZsTa6U2+vQwCg1u5/wt6n9f7ZLkQSvgvadci3rle3J2MfSJi8uJo2myZzl4wbNSDmdYbb6kwzRwAo6sCZ9ye7DUYOHcm3/rz9syPq+zz7vLb+nwijrBSfE7uh9CNfsfOn0AtPUmf96fx5kOy/y/ImtHMMToj1BBLZXUANstbGuObjNJm7ZGbDajGvMyUrvD/h+2d/AsDtS4LPX1EnwKXVeb88/LX+nf1QhVda+3+Oz9bob89Oyydzr89Ik9f+/UikgYUkdkPt4rOToeZRd7t9vHDpobmkmgxuK2GdYC/pwZDEbJrMXTBodD+mnRB8LuX2u352tM6zDi2k5vqdYe3Tt9dItrRtypCHn6SICfxHWO9A6BefnLt/TlgxeOnqU9MZ2uVJr8NwXM4Xl4dzsN/pc0Wnv19F/n5FNj28yNWY3PDmFfew8aJ0JEYjiILRZO6CF5v9u8Btyhw+TGnj3HC58+Ys5NH+kc+TftU2Z75cPrn6Xs7xmRFPgGJFA98er9bayC8TB/LNhfIvn0vH3U5lTn98xaHkWySqq2S9SUrlK1bkm8tv96TuZKHJ3CNpEVwMEky0aSDFvnz5wj0zc9ZVz9xolR3g8x1qX3iPASPyfXFtaduULW2b0q/XiJBjDKW+sgf8X0GZPTzOyVQVLxP7ljh6FLDGST82812KhDAE0mvnr0y+LhevaTJ32LtvTSxwm6u3/UD97QccHSbXZ4Cbo1X8p8CQ4rc3KW0C32Th5DAuPiqoTt+E7y9qp1q/kXGn1XtT1UbcuOE7OlY5jceeCO/XWYd8s1Y6f3zqZKWRajJovXxZzrpu3aMbIHBeQaNuHOz2uGT3r5x+ZIlj5bnFiTsNKR9jq1UpcJtXbvoPAOOmhT56omLWDnamVAKgiDlGhhSNLEAHhZTM7c9UsAtHLlz+Jysan+xQVP45eWHNmYf/BNJDKjEWnRYXXHQ5F3B5zuNQ0/FN66cw9o6e7gTl4/6HenE/wEXpvBTR3Zryq7pnD5R2pKgCTbr2IavOILEXT7VSaQkP58PRlrnD1hapk2/dHSu/5JEf/i/f+nASzGkHjl86XTvCeUbKmD0FbuOboGN1ZWGHJu3yDQvcvzuy26UFHDHi4Eu54e/Qp5gNp51bNNPZ7pHIXvLxvapmBe8KiZdupnhwX/Or6chCxrbw7l72mswdFKiL5YT12+gzaHxIZQSclMmBZNRp6scBnws2DttfgiwbwheDP+fun5OrXx4g/fzzqB9BH2rjo0tD3rbKIWsagDO2bAi7nrw6dXnc7/qC3qKCkl+5vQcjisctfTavC3ufhtutkU8nRPhlHCrvx47kVrRIUca0vZNqJ5zoWQyazB20dF/+uSO6/fgGfQePjbjMm9ZPAUim+sQAABkSSURBVJxpBZUoWarAbSTXsglYd4Mja8LqM8921p8reP+aB/JtVuxo+ONzL1oQ/Ce7b9XlDx6g37zPeO4u97sVolEt4/iJZ6eFO2z0xtvuD7uOvu3vYOAfX9Ln8acK3DbYcNhIZU/+Znze/VYH57ky/3wwbYuH/0UYLU3mDprSMP8twzKOHc1ZTstzqW+kJ0Bj3yoJUGNI4VsbnXTQahGnHvF/ufPjA57J9fhfN97OdZunhhpgSB5+dKCj5TlNjOHB7oOY1+5KV8qvdCD3iKIKu5yfSfLEqlV5sFvfkLYdsW4WQwq4YC0cYrL8fqY+u/LugPPPu9XD/WqLq10qOTA9AeqgtUXq5luX4vPH9dCM9zhWshREcJ9I3/7rrBin80A5OzuK8/fNAvy/ptq7rPlBLljyN83NKnoPDP1Xygu3PMonPiedBMNZG/9mab1Tc2137v45/Fo6PWhZFfYcyLeuTKnKIceSWAKNPspt4H+G+NkojL+tKEeM3HbXwyFtVyNzA7f+MSeiz02onP5EFSuS5nCJBdNk7pBub42CWu3zrfe9fiNvIgvlBKO/sdWhJvNb13zDO3VDOyET7m+E1KysnH72OtsDz7lRzm79des1LOA2oTIIt556IW/nmQLmk6vv5cJvP+CvtFMC7jv8gf5R159XsF9Wpx5dli+eWM0zfrzC2FbnljqHttCjx5NehxH3ou5mEZFaIjJdRJaKyJ8i0s2JwBLNB34SOUC/p16MqlzfZH7/ok9JP7SAIiH+OHw2gtn7/F+YkzsJXbz7V1rNis19Djstn0xJc7xV3ezc4NMkxI14O0PnAre/K8rttd73Gn5mf/RX9xlHlwR9Ptk50WeeATxqjDkNOBvoIiKNHCg36WW37C7eHXxCKLBGlAx6ZBCTr7jLnW6WYK3GPE+9c+1D9Bs8JuAFOC3mLY4qFN/ZDp+6vx9Xb/jF73aJcPf0cCL08vWkH1pgxxA/aXDgI4N4dN4nDL7yHr/P37/kc3ou+jzncamMIzGZrz5eRZ3MjTGbjTHz7OV9wFKgRrTlJpLh/btGtX/jBcGuLsv/x1n7YOjjnMN3vL5go1kAzpm7iFYH51Flw45c6/sMeO74TIsRfLa6/DjJ7/pQPqhOdWUkwheFkyZfcVe+eW4KEosj9PijgylfoYLf5wZ1GUCPR3yvJC28iRwc7jMXkbpAM2CWn+c6A50Bateu7WS1nttbq2pE+4WSME5ctxXqQNPla3LWtVy4itNLrnX0lnSRfAxyphBwePBF30Fjwt7H3xdPof1oF4LvoULwEsPmWDIXkdLAf4H/GGPy3dbaGDMRmAiQnp6eVJ+zN052Mpvl/jN9YsBonoBcZ/IfHWDNvzEhikuj62SsJdAIlLzitpUauwnCAXhw/sfsrlCGUI+b16LpcrButhD8dXr5IS56LMPPWvG7mO3J1T/yZP0Lcq0zee4BELd/6yFwJJmLSFGsRP6OMeYTJ8pMZP/aOp3PqwS4U46PWPTv3brmG0ocPJJvWFezHatyPa6yYw9UgZo7/oFyroeVkAb2GBr2PpG+w91+nkQqKa4Ox/MvdzJ7+PcP2F25LFPqNI9xHMGVPRDg9mxBDvgD93Rj/NTvc+Y48pX3s5iakniX4ESdzMWalf3/gKXGGOen7otzw4d0h/PuzHnc6uA8JnbszucOTSgUrUAjWvK2PwY//CTFB3enz4DneNOOPXubSBrARU12y8m5yzLyhZGS/5ObGubdlqJx6pbNzKkXINn6CaPm7l3MCjy1ey5P9H828sByQsh9xCJpc/brZV3J2WRq+Bf33P3XF+wvXZxE+SWT6Jz4+mkN3A5cJCIL7H9XOFBuQvihxVm5Hn925d2OlV0uhvdMBOvEJcDNa6fQKNe8J+GngStm/cJlO2dwefNLHYrODz8Js/jh/Pchdcszdz8R+Ek/h6xb08vp+vt77gWULU46MYc/2J9xt/ufx8Y7cXJwXBB1y9wY8zOF+HzEH8VOz1n+96bvCacVUlD/XM2jod+uzUnZ85e0+P4LIMKTo/2fcyyeUPoxi9jDGcXkv2CnwbGVeNU69I28YePG9G08kvF5frU1PLaC5UULngK4/a6fmVLhPMCa6GyvnFBwpXnUyNzAxtSaBdYViuRNi5CZ5d1UtpFKvI6hOFbt78S5e0poH8Q4+7jmCcf3vottZ82l/T+/0K5Ri+P9n/bzKR7OMR2KtCxrwqmCzqG8dX1ol79n81de723Bb9OXyCcAff9AIr0HiVvnsWa0PLXgjaKkl/M7qG+/8Po5E+UCh5hfhh6BPgN9Ttd8+36u5+I/+hgxhhtuuYeuIZzPCeWOTF4e10C3EIzXz9TJpYq7Xoe2zFXCis+Pbfgi+bKMh9ceDzGo4zSZR2Hw6OOTN3X+M8BNJYLw95PWpMRPO/LsLcsBKH/Y45smhHlIErurIHoxS7IRfAkNWvItw5bFdm7xYLIkvrvgwqHdLFF4odn1OcttqjfzMBJ3jL/tMcYDD753/MOXZmI3WiQ/d9JUcXOIw1LClbLjSdn9hwrcxu0etfu79HK3gmxhv5DE/52hLXOHXHJdh7D3Kah/L976qk/KWM0DP79f8IYxEuj4hNtvet/v/+PG9d/lWtdxjrM3xvBa58Wf8nC7W7wOw3PtNlnnC1KNvytI3fX66XWZ2Liua+VryzxC3/33M6hg3by5xcH5FIYLI+oc2JozFj2W3G4z9e09Mme5fNZO2m2eT99eI4PsUbAaB3byVzlIC3Sloo9wXp+YTIyk0uTAcn4p3SLk/QZ3HRRGLaGpuWEHeHfLy4iMvf1xxgJ1p/1KZozT3+WV3b20WlvmEfro2Oqc5Vp7YntxT6xlj2wIZYSDq3EUuCLQfrnjvnfJ/2h2+A+Mn9ez9OKLGX/bY5EF6OP6IyV4aP7H9Os1It9z9y75jPsXfZo/zhCuXi2LdbOPFguXRR1jNO5a8SUDHg1/egPlHm2ZR8h37pXqf613rNxwTt5dvnMGX1ds41jdgdRZuo7zS8+i0fxlcK3r1bluaJcnXa/juo6duC5g/da9SF+OYsqHQEPzIDYngEd0Du0+n64p4HuvoCko4qsD0xmazB3Qb3D4U7aGourB3UGff/3fXRk6qjcpR/JPpFWgMPrjew2xb3f3r/CqiBeFdXRLsIRf4L4hdv7UzFzPUYn9/S4DcXOceRmzh0aHVjGrZHiTjtVL2UwsumE1mUep4bHluPFGFTeHOGvphgJbwv165v8ZH0iRosWijMpdt6z5hjL7DuX7Yqqwax/UhFq78t8+THlrTrvY34XeSWkpVgoslVXw8NsVF10AXEDVMH9R/XbB5ZGEFjZN5hEYPLo/2MMSr545G9o7d5OI0tv+gVPhtqVT6NEn+psgJ5LRAWZ4HNRtMMWHPs5lbQN1XCgVmZPrNeKmdd/TvkRlr0OJmibzCHx7xvGW4+MDno64HH8///sMHE0f8GAe6/j2RL/8x7nI4WMeRJL4ypg97As0SVchNPbW6E94xwMdzRKBVUVO8joEhf8E70/s5utIjAtP7pv5JQ8u+NjrMCJyUsaqoM+3Wb0CgKZVasUinLiiLXMPxeukQCq2IjlBG8oFZYFOgPbsMyrs+uJFzvDYAC9/0IMDsEbUF75fttoyD9PQvg/mLHed/UHA7a7ZOj3fVYUq+dTJWONgafrlHo1wRu+cfNS6TuSEcuVdiib2tGUeplcvvj1nuW/PpwJu91LH7rEIp1Bq/88vnLRuS6E+r+DWBVwJ83USZaBjazTih5nfc1rb5OgvB03mYTssId7EMc4lzIfWj7eu6+J1CDlabV7B2lp1vQ7DccneBXjaGc057Yz4ukl1tLSbRSWtWKSj4kfib0RNJDfgVolPk3mEct/wWMWT9PVWf2jlPdY8JonSyHQ+ByfIC1eO0GQehikfH78BRYe5iz2MRAXz9L192NK2KaVCmLEwaknWCm54Yk3KmL1cumK+J/XfvewL7vnrc0/qTnTaZx6Gr/cvg4p1AejmwNWZYt9oWCSVouYox2I4x0WS5aCEFk9zx9x4VxduBLjofE/qH/5A/4I3CuDsNYlzQ3U3OJLMReQ14CpgmzHmdCfKjEfv1bnU0fLu+P1bljesSctm7cj85UMOVSobsxEa+gM8DumbUqBg4+t7do38iyAZONUyfwN4HnjLofIKhX69jw9tvLjDVR5GohJF3nxfPXMj+S+QiZ+WvluS/xWGz5E+c2PMT0By36HBR6flk70OQSkemvcht82e4XUYKoBixOCcjQ/tMw/RiH7d4eI7ASizaavH0aiQxKL55tNUrrk5wvZMhHEOeHR4ZDuqmJjStApzd66NWX0xS+Yi0hnoDFC7du1YVeuY5U3r5Sz3GTTew0hUuGJ1AcygR6K7z6Z2HSSXU8rX45Ty9Qre0CExG5pojJlojEk3xqRXrpx4cwfPqFB4Lx1XwehZSxUfdJx5iPZJWQAu3DPT40hUfIm/9nSrhUtpdHQpp+yOv6tT4901m35FTBbVayZe74EjyVxE3gN+A04RkQ0ico8T5cajxrP/8DoElWQa7dgAQPHdBxwpr3ffZ5h26c088J9+jpRXmIy57VE2X9Sc8hUqeR1K2BzpMzfG3OxEOYmg/7AXvA5BJZkeza6g3n9f4/EQb7ahlD/azRKCrpP0Q5aIEqU3u8Fpp2kij0LlrG2UMvu8DsNzOjQxBD9XbwxAxaztHkeiVPz10Xtt0cXtvQ4hLmjLPASbU6oD0HH+dI8jUUop/zSZh6FlnRZeh6CUUn5pMi+A7z0/L73hWg8jcVai9CcrpUKjfeYF+P78C70OQamI3bL2G1KyTNLcL9Wte58mA03mBfgr7RQALtg7i/yz0yUuPY1WOIy+q7fXITiq4/LlzKi1k07XdfI6lLijyTxETWb/AR28jkKpwu3hh/vwsNdBxCntMw9i0usv5iz3HTrBw0ickVaieM6y/lhVKrloMg/i6dr1vQ5BqVzO3WVNJ1E5AS83V+7SZB7EtpQqADQ6utTjSJSyDDzzMt47tpG7uyRXX3gi6bDnV5458pfXYeSjfeYhmHZpoZl6RsW5eiefTL2TT455ven/LOOrilViXm88evmah7wOwS9tmQcwbGQvr0NQUYvBmJ1CMixo2DnXM3btr16HoYLQZB7A+JbJ3RqvvW6b1yHETKzuNJTMqtWowU13xWeLVFk0mReg69RXvA7BcWIy6dfzKa/DUCG6bOcMHpr3sddhqDinfeZ+DB/cHdpYN29OhiGJKrG98e+uXoegEoC2zP0YZyfyZFVIunmVKlQ0mQfR5buXvQ5BxTvtjldxQpN5Hp3ffy5nuf/wF4NsqZRS8UOTuY+MjAw+r9IWgId/f9/jaFRC0D4rFSf0BKiPpj9Og5QTAejXa4TH0Tjv6KHDXoegfHT9/T0yixRNmulplbc0mdsGjhvIjjOsm088MvPdpP6A6bjr+NC310ivQ1BJRLtZgEGj+/HyGcfvItTniVEeRqOUUuFzJJmLyGUiskxEVopIQs0A9OB7o3mx2b9zHm9J4ha5coH+yFFxIupuFhFJBSYAlwAbgNki8rkxZkm0ZbvhjVcnsHHbWn5udjrzizeBqhcBUMrsY9VFbTyOTimlIuNEn3lLYKUxZjWAiLyPdU8ex5P5wHED2V6pLFkpKWSK5Pk/hcwUsf6XFDJTUshCyExJIZNUDqWmsT+lFDvqn8XRk1oDVgKvkbGZ9vPm0K+3dq0kG5OiQ01U4eFEMq8BrPd5vAFolXcjEekMdAaoXbt2RBUtrVudn8rkLlpMFqlk2v8yji+bTFLJylkunnWEE4/t4tQj66i8dx8Vt++lRfVmXH7zjdD+xojiSTQZx454HYInNKWrwsCJZO7vs5KvJ9EYMxGYCJCenh5RT2PjX+fTtMhijmUcAzGULFmOrt0fIzU1jZSUVFLSSiGiH12lVOHjRDLfANTyeVwT2ORAufkMHPGSG8UqpVTCc2I0y2zgZBGpJyJpQEfgcwfKVUopFaKoW+bGmAwReRj4FkgFXjPG/Bl1ZEoppULmyBWgxpivgK+cKEu5JzMzw+sQYirt8FEAyh/Z73EkSrlPrwAthArL5fwDuw/lzpVf0mbJ316HopTrdG4WldRG3tfX6xCUigltmSsVDR0Jq+KEJnOllEoCmsyVUioJaDIvRM5ueSmnHlvGncu+9joUpZTD9ARoIXLeZZfwg9dBKKVcoS1zpZRKAprMlYqGPWS/yZHF3sahCj1N5kpFwx6aWCQr09s4VKGnyVwpB+hwc+U1TeZKKZUENJkrpVQS0GSulFJJQJO5UkolAU3mSimVBDSZK6VUEtBkrlQUyu7YB8Ap2zZ7HIkq7DSZKxWFfr2foseMdxh9V2+vQ1GFnCZzpaLUc8DTXoeglCZzpZRKBprMlVIqCWgyV0qpJBBVMheRG0TkTxHJEpF0p4JSSsWnu5d9Qaflk70OQ/kR7Z2GFgPXAS87EItSKs4Nf6C/1yGoAKJK5saYpQAiOgGoUkp5KWZ95iLSWUTmiMic7du3x6papZQqFApsmYvI90BVP0/1NcZ8FmpFxpiJwESA9PR0E3KESimlClRgMjfGtItFIEoppSKnQxOVUioJRDs08VoR2QCcA3wpIt86E5ZSSqlwRDua5VPgU4diUUopFSHtZlFKqSQgxsR+YImIbAfWRrh7JWCHg+E4ReMKj8YVHo0rPPEaF0QXWx1jTGV/T3iSzKMhInOMMXE3dYDGFR6NKzwaV3jiNS5wLzbtZlFKqSSgyVwppZJAIibziV4HEIDGFR6NKzwaV3jiNS5wKbaE6zNXSimVXyK2zJVSSuWhyVwppZJAQiVzEblMRJaJyEoR6e1yXbVEZLqILLXvptTNXv+kiGwUkQX2vyt89nnCjm2ZiFzqZtwiskZEFtkxzLHXVRCR70Rkhf1/eXu9iMg4u/4/RKS5Tzl32tuvEJE7o4jnFJ9jskBE9orIf7w6XiLymohsE5HFPuscOz4icpZ9/Ffa+4Y0qX+AuJ4Wkb/suj8VkXL2+roicsjn2L1UUP2BXmOEcTn23olIPRGZZcf1gYikRRHXBz4xrRGRBR4cr0D5wbu/MWNMQvwDUoFVQH0gDVgINHKxvmpAc3u5DLAcaAQ8CTzmZ/tGdkzFgHp2rKluxQ2sASrlWTcK6G0v9wZG2stXAF8DApwNzLLXVwBW2/+Xt5fLO/RebQHqeHW8gPOB5sBiN44P8DvWnERi73t5FHG1B4rYyyN94qrru12ecvzWH+g1RhiXY+8d8CHQ0V5+CXgw0rjyPP8sMMCD4xUoP3j2N5ZILfOWwEpjzGpjzFHgfaCDW5UZYzYbY+bZy/uApUCNILt0AN43xhwxxvwNrLRjjmXcHYA37eU3gWt81r9lLDOBciJSDbgU+M4Ys8sY8w/wHXCZA3FcDKwyxgS7ytfV42WM+QnY5afOqI+P/VxZY8xvxvrUveVTVthxGWOmGGMy7IczgZrByiig/kCvMey4ggjrvbNblBcBHzsZl13ujcB7wcpw6XgFyg+e/Y0lUjKvAaz3ebyB4MnVMSJSF2gGzLJXPWz/VHrN52dZoPjcitsAU0Rkroh0ttdVMcZsBuuPDTjRo9g6kvsDFg/HC5w7PjXsZTdivBurFZatnojMF5EfRaSNT7yB6g/0GiPlxHtXEdjt84Xl1PFqA2w1xqzwWRfz45UnP3j2N5ZIydxff5Hr4ypFpDTwX+A/xpi9wIvASUBTYDPWz7xg8bkVd2tjTHPgcqCLiJwfZNuYxWb3hf4L+MheFS/HK5hwY3ElRhHpC2QA79irNgO1jTHNgB7AuyJS1q36/XDqvXMr3pvJ3WiI+fHykx8CbhogBseOWSIl8w1ALZ/HNYFNblYoIkWx3qh3jDGfABhjthpjMo0xWcArWD8tg8XnStzGmE32/9uwpiFuCWy1f55l/7Tc5kFslwPzjDFb7fji4njZnDo+G8jdFRJ1jPaJr6uAW+2f1djdGDvt5blY/dENC6g/0GsMm4Pv3Q6sboUiedZHzC7rOuADn3hjerz85Ycg5bn/NxZKZ388/MOae3011gmX7JMrjV2sT7D6qcbkWV/NZ7k7Vt8hQGNynxRajXVCyPG4gVJAGZ/lX7H6up8m98mXUfbyleQ++fK7OX7y5W+sEy/l7eUKUcb2PtApHo4XeU6IOXl8gNn2ttknp66IIq7LgCVA5TzbVQZS7eX6wMaC6g/0GiOMy7H3DuuXmu8J0IcijcvnmP3o1fEicH7w7G/MlUTo1j+sM8LLsb5x+7pc13lYP2v+ABbY/64AJgGL7PWf5/mD72vHtgyfM89Ox23/oS60//2ZXSZW3+RUYIX9f/YfhQAT7PoXAek+Zd2NdQJrJT5JOMK4SgI7gRN81nlyvLB+fm8GjmG1cu5x8vgA6cBie5/nsa+mjjCulVj9ptl/Zy/Z215vv78LgXnA1QXVH+g1RhiXY++d/Tf7u/1aPwKKRRqXvf4N4IE828byeAXKD579jenl/EoplQQSqc9cKaVUAJrMlVIqCWgyV0qpJKDJXCmlkoAmc6WUSgKazJVSKgloMldKqSTw/2IE2oRL0J8FAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 100     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                mine_list[i].step()\n",
    "                dXY_list[i,-1] = mine_list[i].forward()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current results saved.\n"
     ]
    }
   ],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    mine_state_list = [mine_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'mine_state_list': mine_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current results saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth is 4.982193620464953 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = g.ground_truth * d\n",
    "print('Ground truth is {} nats.'.format(mi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "mi_list = dXY_list.copy()    # see also the estimate() member function of MINE\n",
    "for i in range(1,dXY_list.shape[1]):\n",
    "    mi_list[:,i] = (1-mi_ma_rate) * mi_list[:,i-1] + mi_ma_rate * mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gUVdvA4d9JgRBKAoTeEor0HqSqQUEwNEUU/QDFAr4iKioWRGz42iv6WrCBgBVRRJqCgtIJEIoUCRAlgQABEgjpyfn+mMmym90kC+xmFvLc17UXszOzM89Oln32zGlKa40QQghRmJ/VAQghhPBNkiCEEEK4JAlCCCGES5IghBBCuCQJQgghhEsBVgdgLywsTIeHh1sdhhBCXDQ2bdqUrLWu4Y1j+1SCCA8PJyYmxuowhChVB1MPAtAgpIHFkYiLkVLqH28d26cShBBl0agfRgGwYvQKawMRohBJEEJY7Kkrn7I6BCFckgQhhMX6NO5jdQhCuCStmISw2P6T+9l/cr/VYQjhREoQQljszvl3AlIHIXyPJAghLPZc1HNWhyCES5IghLDYVeFXWR2CEC5JHYQQFtuTvIc9yXusDkMIJ1KCEMJi9/x8DyB1EML3SIIQwmIvXvOi1SEI4ZIkCCEs1qNBD6tDEMIlqYMQwmI7ju5gx9EdVochhBMpQQhhsfGLxgNSByF8jyQIISz2Wt/XrA5BCJe8miCUUvHAaSAPyNVaR3rzfEJcjLrU62J1CEK4VBoliN5a6+RSOI8QF6XYpFgAOtTuYHEkQjiSW0xCWGzCkgmA1EEI3+PtBKGBX5RSGvhIaz298A5KqbHAWIB25ctDVJTjDjffDOPGQXo6REc7n2H0aOORnAzDhjlvv/deGD4cDh6EUaOctz/yCAwaBHv2wD33OG9/6ino0wdiY2HCBOftL74IPXrAmjXw5JPO299+Gzp0gGXL4IUXnLd/9BE0bw4LFsAbbzhvnzULGjSAb76BDz5w3j53LoSFwYwZxqOwRYsgOBjefx++/dZ5+4oVxr+vvw4//+y4rUIFWLzYWJ46FZYvd9xevTp8/72xPGkSrF3ruL1+fZg921ieMMG4hvYuuwymmx+JsWPh778dt3foYFw/gJEjISHBcXv37vDSS8byjTfC8eOO26+5BqZMMZavuw4yMhy3DxwIEycay4U/d1Bqn723W0wwPhszCsUgnz1jWT57ztvtP3te5O0E0VNrfUgpVRP4VSm1W2v9h/0OZtKYDhBZubL2cjxC+JwOVVtCWiWrwxDCidK6dL6TlVLPAmla69eL2icyMlLLnNSirNmYuBGQympxfpRSm7zVAMhrHeWUUhWVUpULloFrAekNJEQhj/76KI/++qjVYQjhxJu3mGoBPyilCs7zpdZ6iRfPJ8RF6b3o96wOQQiXvJYgtNb7gfbeOr4Ql4o2NdtYHYIQLslYTEJYbM3BNaw5uMbqMIRwIv0ghLDYk8uNJqrSD0L4GkkQQljso4EfWR2CEC5JghDCYs3DmlsdghAuSR2EEBZbGb+SlfErrQ5DCCdSghDCYs+seAaQOgjheyRBCGGxz4Z8ZnUIQrgkCUIIizWu2tjqEIRwSeoghLDYsv3LWLZ/mdVhCOFEShBCWOyFP4yhuPs07mNxJEI4kgQhhMVm3TDL6hCEcEkShBAWaxDSwOoQhHBJ6iCEsNiSuCUsiZOBjoXvkRKEEBZ7edXLAPRv2t/iSIRwJAlCCIt9Pexrq0MQwiVJEEJYrHal2laHIIRLUgchhMUW7FnAgj0LrA5DCCdSghDCYm+sfQOAQc0HWRyJEI4kQQhhsbk3z7U6BCFckgQhhMXCgsOsDkEIl6QOQgiLzds1j3m75lkdhvBxuXn5pX5OKUEIYbFp66cBMLTlUIsjEYXdO/1P9h85xdIpA1xu3xp/HK01HSLcLwUujT1I9cpBRDap4bRt8/5kJs1ZD0CLeqG8c2dPhzi+fqgPVSuVP493cn4kQQhhsfm3zLc6BFGE/UdOAbBp/zE6N3b8Qr9t2m8cSc1wWLf4qWh+XH+AHzbEk56VQ1pmLt8/ei2VggIBmDhzLdv/PQHAj4/3o0K5s1/BWmtbcgDYnZjiFMctby0rMll5g9xiEsJiIUEhhASFWB3GRSMvX5ORnVvifjkXeEtm4sy1tuUn52ywLa/ZncSqXYedkgPArW8t46Nfd3E0NYO0TCPGG1/7xZYUCv4F+L+3l5OXb8Q4a+Xf/LY90el4efn5bNp/zGHdgSOnOJGW6XJ/T5MShBAW+2bHNwAMbzPc4kh8y6fLdxOXlMpLI7o6rB/00mLy8jWfjYuiXvWKLNuWwLZ/jjNhYDv8lAJg075jPPnlBt4f04smtY3k+/bP26hWKYjboi4r8dynMrIdvszBuM3z1ujuPPfdpiJfl3Im2+X6xONpNKlVxWFdelYu0f9dzJRhnZj9x17b+qoVy3PyTBYA0f9d7HSso6cy+M/0PwGIqFm5xPdyIaQEIYTFPoj5gA9iPrA6DMs9/90mbnhlKQCxB5L5ds0+Nu9Ppt/UhRw/nWnbLy9fA7AlPpnVu5N4bf5WlsYmsHp3km37k18av/jHfbyKHf+e4ONlu1i85SBz/tyL1rrYOGIPJPPw52tsz69oafR033/kFI/OWuewb68WxrYRVzRzeaxqZn3BWz9v54ZXjff27M2RDvtMnbvZ4fkXD/Rm0tCOTsfqaNZzZGTn2dadysgp9r1cKClBCGGxRSMWWR2C27TW9H/BiPerh66hWqUgl/sdTc1g8pcbeGpYJxrVKPlX7sa4o7Yv+H5TFzpt/7+3l7N0ygCHbUkn03l30Q7b8xfmbmZY98bMXbvf4bWP2N0qAkhNzya04tmKXq01q3Yn0bVZTcoF+PP47LP1AD890R9/P8Wfu4xf8n8fSnU41lPDOnE0NYNaocHM+fNsKaB9eHUejG5LveoVnd5Px4jq9G1Xn1+3Jbi8FuUC/AmtWM5hXbM6ITwyuB0j3/nNdp0AHiuUsDxNShBCWCw4MJjgwOBSPWfC8bTzet3OhJO25b8Onixyv1HTfuPf5DTGfvhHicdcvOVfnvpqY4n7zVt/wOH5d4USAeCUHFxJPHHGtrzlQDL9X1jEC3M3M39jPCfTshz2LR/oT4C/H73b1HV5LKUUtUKNv91d17SwrX91VDfqVa8IwIMD2trW169ekaByAUwc0p537+rpdLzvJvYFoFL5QNu6iYPb897dvahorttR6NaXN0kJQgiLzd42G4CR7UaWyvkKftE+NLAt/Ts2PKfXrt1zxLack5tPelYuiSfOUL96RW587Rd6tajNFS3rFHuMrfHHaV43hCCzBc/bP293uV/jWlVsrXcAPvplp8v92jWqxrFTmRw+mW5b17NFbZ6+qbPL0sjhk+m0blCNvHzNE3alhdMZOcTsO+a0P8ATN3Rk/5FT/HMsjavb1KVGlQqMvMrxttLNPZpwedOa1KnqmOyjOzXkw6V/kZWbb2u2CnBZ3VDKBfiRnZvPosnR+Psp27aCVk8ATWobdRdB5fxR4FYFvadIghDCYp9s/gQonQTxzeo42/JbP2/nilZ1bL9MC/v7UAo5efm0blANgA17jzr8ak9Jz7bdVy+wcudhVu48bHveIby6w/bt/xy33RYpfMvovbt70SCsEvFHT/PgZ6t5f0wv5q7dT2z8cYcv7u8e6ctNb/xqe346I4eXR3Zl9Lu/o4HP74uibjXj1/v7Y3ox7uNVACyY1J9BLy1hd2IKfdrVJzY+udC12cctPZsAcH90GwZ0ckye/xwzSl2tGlRlUGS4y2sWXkSl8awHr2FXwkmHL36Abx/pS3pWrkNyAKgSfPYWU22zhOKnFBXKB5CelYufUuSXUJfiCZIghLDYr6N+LXknD5j63SZW2d2/Bhj66i8sfiqa3Lx8ygX4O2y7/9PVgPFFnp2bx5SvjdtAnRqHsTX+OD9ucLzlU1jL+qFs++cET85ZzwMD2lI7NJhldk0zM+1+CVco50+zOkZroxb1Qm1t/W/q0YQbuzfmuhfO1tNUCS7Hf65txYe/7OSVkV1p07AaAf5+LHHRPyC8ZhX8lGJM35a291fQPNW+6WqBr1fvA2BAp4YopZy2A3RtVqvY9+1KSHA5ul3m/LoK5QIc+kIUCC4fwDt39iSiZmXKB579u1Q0E0TNkCA+H98bP6VQT59zOG7zeoJQSvkDMUCi1nqgt88nxMUm0N/1L/hzkZSSTvkAf/LyNWFVXFccF04OBey/fAu+mFftOuywz2NfnK0MfWlEV259axlHUpz7AdjblWB09Nq032htdGO3xizZchAwksdGu1LBj48XPZuen1I8Mrgdb/y0zZZEbugawQ1dI4o9P4C/n2LxU9G25/WqVSQzO9flrSd7rpLD/Mf7kZ6dW2TFvKe1qBfqtK5yhXIcO5VJ9cpBtia93lQaJYgHgV1AlZJ2FKIsmhE7A4DRHUa7tX9SSjrb/zlB3/b1AaPN/u3v/u6039IpA1i/9whLtxxktV3dwfjrWpObp9l35BS/bnVuSXMmK8eh6WXsgWR2mb16o83bLsHlAjiBUaH77l09WRDzD+Ova8Pgl425tf93dy/u+2SV7RiB/n6ctmuSuSshhYPJxi2bLydcU+J7vrZ9A5rUCqFxrQtr95944oxDJTUY1+nNBUZT2eIElQuw1ZtYpUI5ozRRvXLpJCmvvlulVH1gAPBf4GFvnkuIi9X0mM/IyslzO0EUJIOeLWpz8kwWd/5vhcv9TqRl8vTXMQ7rhvdo4nD//N9jaew5dHZIB1e/rO2bfRa0yEmw+5K9rG4ojww2fu1+cu9VVK4QSGjF8rw/5grGfWx06MrKzXNqOTVzxd/A2b4CJSmorPWkbx7uA8BDA9uRm6dZvj2R6y8P9/h5PKWg5VhpJQhvN3N9G3gMKLLPu1JqrFIqRikVc+yY6xYEQviyDXuPkpmT57Cu39SF3Gv2di3sTFYOKWeyOJicxqQ566l8+HHCkp+k39SFHElJp9/UhfSbupDE42ecXvvXwbNNHG94dWmRyQHg1reWO627vXdzh+f2yaGwp4Z1KnJbgFmp+txwx05fDcIq2foYNKldhUWTjds7mdl5/H3Y6ENQuOK6qHv93jbnwWtssSqlmDCwLS/c2oV7+7W2JB53dG5sdJbLKvR58xavJQil1EDgqNa66H7pgNZ6utY6UmsdWaOG8+iGQviyflMXMuXrjdzwylJy8/JZszuJmb/vAYyet6npjkMvaK0Z+uovDH9zGeM/WcXm/Y4taW6zu1V05/sr+GTZLvpNXcjf5hf5wzMcO30VGNC5IZc3q1lknFe2qsPSKQOcWssUmPOg422eiJqVaVC9ksO612/rZlt+9bZuDOzckK7FnBOMOoDygf5kZOfy/pK/AKMTWYFr2tYr9vWeNumGsz2UC9fVlAvwp0vT4t+P1f4xb8st2vxvqZzPm7eYegKDlVLRQBBQRSk1W2tdOo29hfAy+1JDvtYMeNF53Jyb3/jVVvGbciaL4W8uc3p9gjbu29dXzhW1Bc1K7/90dZHDOdzdpwU3dW9iiyMvXzPQLpbiRv/8bmJftv9zwunLsnebeg5NNls3qErbRtXtnlezNX8tSYVy/g7DQ9if66rWxfeZ8LReLWtzR0pzrnejgtsXTbuzJ//39nJmju9dKufzWoLQWk8CJgEopaKAiZIcxMVid2IK363Zx1PDOqGU4kxWDsu2JtCoZmXCKgdRPtCfke/85taxtNa8Nn8ry4sYfTNJ/0mbBtUg8WyCuLVXU75aFeewn/1QDkueiubuD1aScPyMQ8c0P6Xw8z9bShjarfgvwioVytHTHE/I3g1dw4Hik4u7KpQLcOjcVdVumItGYd4dbK6wAH8/bunVtFTP6UnVKweV6nDf0g9CiEKOn87kwc+MPgCPzVrHtn/OfWiDIV3CqVIhkFl/7OV0Ro5DcnhwQFveWWj0Hv7qoWsoH3AtFYMCbRXEt/RswujezalTNZg3F2xzOvaro7qhlOLTcVFFnj+4XADp2bmMjmpe5D6FPTc8ktd/2sp3j/T1aL1AUKA/v+84ZHveyW5ehZqhFTx2HuF5pZIgtNYrgBWlcS4hXMnKyWPtniNc2bpOse3Hj6ZmMGra2ZJBScnhzdHdqRUSzIh3jArhRZOjSTmTRfXKQaz8y/hStB+JdN5j11KxfKAtQdi3qV80OZo/dh7iqtbGuD89W9Qm9kAyd1zdgme/iWGfOexE+0KVvK7MmXANWmuHTlYl6XZZLeZOvNbt/d1lP2/CLT2b4O+n+ObhPpzOyCmVtvzi/EkJQlzS7EcfBXjphy10bVaT52/p4rBfXr4mMyfXITmUxL6ov+SpaNuv7oImiAXDKhw0WyON7t3cNqzFiCuaEVze+O/3/sb3ARjXZRy925yttK0UFMjjZqVqweQ3xVVE2ys4ti9Izzp7e+mOq40B7UIrlncYUVX4Jt/5FAnhBTtcjDi6fu9R8vI1Yz5YyaShHWlWJ4Tpv+7kxw3xtn2+uL83pzJyGG929ioYjmLd30dp07Cq0/AIrm7J+Jt1AQVNU9s1Olupaz9pzYK/FwBGgihKwa/wqFalW6nrCVe3qctvOw4R6C+DR19sJEGIS1Z2bp7DtJFXtqrDHzsP07R2FT5etovEE2cY/8kqlk4Z4JAcbugaQa3QYGqFGr18C0oE5QL8ufIcvqDz8ozB1DbGGf17mtZ2Pa3o4hHOrZ8KK2j33srNlkO+pKCkMLhLI4sjEedKEoS4aC3Z8i8Hjp522bFpw96jtsHlABZNvg5/Pz9ycmNISkknyRwa2t9POQwpDY7j+l9Ij9WCSV8KhnY4l/qAwt4fcwXfrI5zGkr6YlAw5qjcUrr4SJlPWCY7N++8hyxOOZPFWz9v58cN8bYpKAv8sP6AQ3L44bF++PsZH/Wgcv5k5uTZKo4D/BSTv3Qc1dNTt0LcHbfnnXXv8M66d4rdp0ntKjx5Y9E9m31ZwaitF5IghTUkQQhLHDpxhkEvLXEYSdSVTfuO0W/qQqbaTRT/4GerHTqcfbD0L9vy7sQUPrSbWGZc/9YOFbaVggI5fDLdNuxDdm4+J8xZxOY/0Z8Fk4oeVfRcBdl9Ic5+8Ooi91t+YDnLDzgPi3GpKOgt3aWJjJRwsZFbTKJUncrI5qbXHec/+G7NPm7q0YTkU5nM+XMv9/VvTYC/H1qfnXx+1e4kjqZmcPB4GrsTHccPWhDzD3f3aUlQoD8zV+yxrZ//RH+HL2mAhZv+sS3b9/AN9Pdz2vdC2R+vRpWi2/v/dOtPHj2vr2nbqHqpdu4SniMlCFFqtNZOyQHgF3PI6RHvLGfR5n9tPYYnF5qneNS03xwmefn5yetsy0tjjXkGTpzOom61YBZNvs7lF779/MD2dRfeaBZaMejC53kQwkpSghCl5ma7aSLv69+autUqMvnLDfybnOYwzPSXf8bx5Z9nh5loUS/UqdTw8X+uJNDfz9ZjODPbqM+IP3YawFbnUFj/jg1pXjeU8JqVmfPH2aErPr8vyhNv0cl7d/cirISK7tfXvA7AxB4TvRKDEOdLShCiVPSbupBT5oQxr93WjcFdwol08570A9FtHZ4/NqQ9DWsYY/h8ZY7nn5GdS/xRIzmUNEJoRK0qKKXo3vzsGETe+rXfrE4IVUuY72BtwlrWJrgepVUIK0kJQnjEkZR0bnv3dx4b0p5diSncc20rW2ugM1lnZxKrUzWYdnajgn78nysZ8+EfRR53yVPR5Nq1Uip8Lzso0J9KQYGcycqxzW1wYwkD1BUomJ2stIecLuz7m7+39PxCFEUShDhvWmsWbf6XbpfVss1j8Or8rQBUDgqkd9t6bI1Pts0r8H+9mjpNWNOwRmUWPnkd6Vm5BJXzZ9m2RNs4Rd+ag8YF+isWTOpfZPPTiuUDWLY1kfTsf/D3UzSu5d7MY0opqTwVohiSIITb0rNyueHVpYBxb71gGIppi3Y47fvlqji+LDRcdVFj8Af4+1El2OhUFt2pIf07NiArJ89hOItyAUW3MAouH2AbiiIvX1s2Q9n5ennVywA80esJiyMRwpEkCOGS/SB3Q7tFcE/fVg5DVo+3m5DeXSFmEiiJn1JOYx0Vx76jXKv6Vc85LqvFJsVaHYIQLkmCEC7tP3Latjxv3QHmrTtQ7P4vjric+tUq8sXKvzmZlsWmQlNperqPgb1/zWkYm9cN5a07enjtPN7y9bCvrQ5BCJekFZMgOzePeev224aUBlj795Ei939kcDsAvnukL1Hm3AWdIsKoFRrMo0M6uBxzZ/4TnuuhXFhB3cTFOE6REL5MShBlnNaaQS8ZcyJ/9Osufn7yOv7z4R8kmAPMLZ0ygM9+2803q/cR1boujw5pT4C/H9e2bwDApKEdmTS0o8MxU9KzAejZvBYp6dkOcxt7Q5uG1dhyIJkHott49TzeMnXlVACmXDXF4kiEcFRiglBKXQZ8ANTSWrdRSrUDBmutX/B6dMKrUtOzHTqvAQ6T3Q/s3BCAO69uwR29m7td+XvX1c3ZtO8YY/u2onYp/KqfclMnDp9Iv2h7Lu85vqfknYSwgNIljKaplFoJPAp8pLXuaK7bobX2+M+1yMhIHRMT4+nDiiLY916+65oWfLp8t+15w7BKfHzvVVaEJYQ4B0qpTVrrSG8c251bTMFa6w2Ffj3mFrWz8E1an23+eejEGe743wrbthdu7WLr1dwxIowmtavIXMFCCLcSRLJSqgnmvB9KqWHAYa9GJTzq5R+28PuOQ4AxdaZ9cqgUFECXpsY8xzf3aGJFeGXe078/DcDzvZ+3OBIhHLmTIO4DpgMtlFKJwAFghFejEhckOzePZdsSaVYnhMa1KtuSA8D9hfovfP9ov9IOTxRy8NRBq0MQwiV3EoTWWvdRSlUE/LTWp5VS7g12I0rd6Ywchr3+CwC1QyvQutAcxnFJxvSaM8f3plZo0XMUiNLz+ZDPrQ5BCJfc6QfxPYDW+ozWuqD31FzvhSTOV25ePm8t2Gp7npSSYev9vGBSfyLM5qb1q1ekdtXgi25ICiFE6SqyBKGUagG0BkKUUkPtNlUBzn8md+Fx+VoXO3Vn3/b1KRdgjHoKRiW18B2Tlk0C4KU+L1kciRCOirvF1BwYCIQCg+zWnwbGeDMo4T6tNY99sc5h3di+LYk7nMpvZt3DOHPmtO3/ngCgT7v6pRukKNbxjONWhyCES0UmCK31fGC+Uqq71lpmM/Ehmdm5aKB8oL9DyeGdO3tQtWJ5aoUG88LcTbb1BdNpPn59B37e9A8PDWpX2iGLYkwfNN3qEIRwyZ1K6i1KqfswbjfZbi1pre/0WlSiWENeMYbcth+Y7vrLw2lR7+xIpgXjE/Vpd3YynKvb1uNqiyfHEUJcPNxJELOA3UA/4HmMJq67vBmUcO2DpX/x44Z42/OHPl8DwCsju9IhIsxh30cGt6dF/aoMjmxUmiGK8zDxF2Mu6tevfd3iSIRw5E4rpqZa6ynAGa31TGAA0LaE1wgPW7jpH4fkYK9wcgBjEp4hXcKlpdJFICMng4ycDKvDEMKJOyWIggmFU5RSbYAkINxrEQkHSSnpHDuV6XLWNoB7rm1VyhEJT/vfgP9ZHYIQLrmTIKYrpaoCU4CfgErA0yW9SCkVBPwBlDfPM1dr/cwFxFqmHDuVwch3fnNav3TKABbExPPe4r8A6NasZmmHJoQoI0pMEFrrT8zFlUDjczh2FnC11jpNKRUIrFJKLdZaryvphWXdir8O8dK8LU7rF02OBiAo8OyfrTSG0xbeNWHJBADe7v+2xZEI4cid+SBCgdswbivZ9tdaP1Dc67Qxjnia+TTQfBQ/trggLz/fITnMvL83k+ds4OHB7fD3M+oT7IfIkFFXhRDe4s4tpkXAOmA7kF/Cvg6UUv7AJqAp8D+t9XoX+4wFxgI0bNjwXA5/SVoam3B2ecoAAD69L8phn3aNqhPdqSHXXx5eipEJb5GSg/BV7kwYtFlr3emCTmKUQn4A7tdau65tRSYMOpicxt0frATgx8f7UaGczAgrhCieNycMcqeZ6yyl1BilVB2lVLWCx7mcRGudAqwAvDdz/SXgE7sZ3SQ5lB33LbyP+xbeZ3UYQjhx51soG3gNmMzZOgRNCRXWSqkaQI7WOkUpVQHoA7xyAbFesvK15r9zN7Pu7yP061Cfhwe1tzokUYoqBMqw68I3uZMgHsboLJd8jseuA8w06yH8gG+11j+fa4BlwYzf97BqdxIAbRtWtzgaUdqkB7XwVe4kiL+A9HM9sNZ6G9DxnCMqQ/46eIKHZziOg3hlqzoWRSOEEI7cSRB5QKxS6neMvg1Ayc1cRfE2xh3lqa82Oqx7ZHA7ygf6WxSRsMrYBWMBGdVV+B53EsSP5kN4SPzR0w7J4ecnr8NPgb+fO20GxKWmegW5rSh8kzs9qWeWRiBlybz1+23LS56KlgH1yjiZSU74quKmHP1Wa32zUmo7LnpAa61l1pnzsOPfEyyNTSCqdV0mDZUqGiGE7yquBPGg+e/A0gjkUpeXr5m7dh+f/bYHgBu7RVgckfAVd8y/A4DPh3xucSRCOCpuytHD5uI4rfXj9tuUUq8Ajzu/ShTlv3M3sXrPEQD6d2zAZXVDLY5I+IoGVRpYHYIQLrlTSd0X52RwnYt1ogh7D6fakkPjWlWYMEDmWxJnPd/7eatDEMKl4uog7gXGAU2UUtvsNlUGVns7sEvJ16viAHj6ps70bFHb4miEEMI9xZUgvgQWAy8BT9itP621PuHVqC4hh0+ms3p3Ev07NpDkIFwaOW8kALOHzrY4EiEcFVcHkQqkKqWeApK01llKqSignVLqC3MAPlGCdxdtRwO3XXWZ1aEIH9W8enOrQxDCJXfqIL4HIpVSTYFPMaYd/RKI9mZgl4LjpzPZtD+ZapXKU71ykNXhCB815aopVocghEvudN3N11rnAkOBt7XWD2EMxCdKsGjzvwC8NKKrxZEIIcS5cydB5CilbsWYdrRgNNZA74V0aTidkcPsP/bSMSKM8JqVrQ5H+LBb5t7CLXNvsToMIZy4c4vpDuA/wH+11geUUhGA1KYVIy9f8+QcY3bVYd2LnTZDCDrU7mB1CEK45GWSEp8AAB89SURBVM5YTDuVUo8DDc3nB4CXvR3YxWxj3FH+PpzKvf1aEdmkhtXhCB/3RK8nSt5JCAuUeItJKTUIiAWWmM87KKV+8nZgF7PfdxwiJLgcAzs3sjoUIYQ4b+7UQTwLXA6kAGitYwEZSKgIaZk5rN97hB7NaxHgL8N3i5Ld+O2N3PjtjVaHIYQTd+ogcrXWqYWGpHYa3VUYPvttNxnZeQyKlNKDcE/3+t2tDkEIl9xJEDuUUv8H+CulmgEPAGu8G9bFKScvn1+3JlCtUnma1A6xOhxxkZjYY6LVIQjhkjv3QO4HWmNMN/olkApM8GZQF6vt/5wgOzefB2UwPiHEJcCdVkzpwGTzIYqxIe4ogf5+dAiXKSSF+wZ/NRiAn26Vth/Ct7hzi0m4acPeo7QPr05QObmswn3XRFxjdQhCuCTfZB6SdDKdxBNnGNxFKqfFuXmw24Ml7ySEBaQdpofE7D8GQOfG0jFOCHFpKG7CoHcppjmr1voBr0R0kdq07xi1QipQv3pFq0MRF5nr5lwHwOIRiy2ORAhHxd1iiim1KC5yKWey2Bp/nCta1aFQfxEhSjToskFWhyCES8VNGDSzNAO5WO349wSPzFwLQKv6VS2ORlyMxnUZZ3UIQrhU3C2mYtvcaa0Hez6ci8+S2IO25cvqSOc4IcSlo7hbTN2Bg8BXwHpA7p24sCvhpG25YY1KFkYiLlZ9vugDwLLbllkciRCOiksQtYG+wK3A/wELga+01n+VRmAXg8QTZ0g4fgaA5nVD8feTRmHi3A1vPdzqEIRwqbg6iDyMIb6XKKXKYySKFUqp57XW75Z0YKVUA+ALjESTD0zXWr/jmbB9w/q9RwH4YOwV1K0mrZfE+RnTeYzVIQjhUrEd5czEMAAjOYQD04B5bh47F3hEa71ZKVUZ2KSU+lVrvfMC4vUpsQeSqVetIo1rVbE6FCGE8LjiKqlnAm2AxcBzWusd53JgrfVh4LC5fFoptQuoB1wSCSIvX7P93xNc1aqO1aGIi1zUjCgAVoxeYWkcQhRWXAliFHAGuAx4wK59vwK01trtn81KqXCgI0Zld+FtY4GxAA0bNnT3kJbbf+QU6Vm5tGskA/OJCzO6w2irQxDCpeLqIDxS46qUqgR8D0zQWp9ycZ7pwHSAyMjIi2Yiou3/HAeQBCEumCQI4au82uxGKRWIkRzmaK3drbu4KGz95wR1qgYTViXI6lDERS4nL4ecvByrwxDCidcShDLuSX0K7NJav+mt81ghX2t2/HuC9lJ6EB7Qd1Zf+s7qa3UYQjjx5nDfPTHqMbYrpWLNdU9qrRd58Zyl4sCRU6Rl5tC2UTWrQxGXgLs73W11CEK45LUEobVexSXa+3rrPycAqX8QnjGy3UirQxDCJen6ex62xh+nTtVgaoZUsDoUcQlIz0knPSfd6jCEcCIJ4hxprdmdeJI2DeT2kvCM6DnRRM+JtjoMIZzIlKPnKCklg5Qz2bSoH2p1KOIScW/kvVaHIIRLkiDOUcHorS3rSYIQnjG8jQzWJ3yT3GI6R7sTUygf6E94zcpWhyIuEamZqaRmplodhhBOpARxjnYlnqR53RAZ2lt4zJCvhwAyFpPwPZIgzkFWTh77kk5xY7fGVociLiEPdH3A6hCEcEkSxDmIS0olL19L/YPwqKEth1odghAuyX2Sc7A7MQVAWjAJj0pOTyY5PdnqMIRwIiWIc7ArIYVaoRWoVkkG6BOeM+zbYYDUQQjfIwniHOxOPElr6SAnPOyR7o9YHYIQLkmCcFPyqUyOncqkhdQ/CA8b1HyQ1SEI4ZLUQbhpV6LZQU7qH4SHJaUlkZSWZHUYQjiREoSbdiemEOjvR+Nabs+0KoRbbpl7CyB1EML3SIJw066EkzStU4VyAf5WhyIuMU/0esLqEIRwSRKEG3Lz8tl7OJUBnRtZHYq4BPVv2t/qEIRwSeog3PBvchrZufk0rxtidSjiEnQw9SAHUw9aHYYQTqQE4Ya4JGMgtaa1JUEIzxv1wyhA6iCE75EE4Ya4w6eoUM6fetUrWh2KuAQ9deVTVocghEuSINwQl5RK41pV8FOX5BTbwmJ9GvexOgQhXJI6iBLk5Wv2JZ2iWR25vSS8Y//J/ew/ud/qMIRwIiWIEiSeOENmTp7UPwivuXP+nYDUQQjfIwmiBHGHCyqopYOc8I7nop6zOgQhXJIEUYK4pFTKBfjRsEYlq0MRl6irwq+yOgQhXJI6iBLEJZ0iomYVmWJUeM2e5D3sSd5jdRhCOJFvvWJorYk7nErTOnJ7SXjPPT/fwz0/32N1GEI4kVtMxUhKyeBMVq5UUAuvevGaF60OQQiXJEEUo6CCWpq4Cm/q0aCH1SEI4ZLcYirG3qRU/P0UjaSCWnjRjqM72HF0h9VhCOFEShDFiEs6RXiNyjLEt/Cq8YvGA9IPQvgeSRBFKKig7nZZTatDEZe41/q+ZnUIQrjktQShlPoMGAgc1Vq38dZ5vCX5dCap6dlSQS28rku9LlaHIIRL3qyDmAFctDOhxB0+BUBTqaAWXhabFEtsUqzVYQjhxGslCK31H0qpcG8d39viklLxU9C4ZmWrQxGXuAlLJgBSByF8j+V1EEqpscBYgIYNG1oczVlxh1OpX70SQeUsv0TiEvd2/7etDkEIlyz/9tNaTwemA0RGRmqLw7GJSzpF+/DqVochyoAOtTtYHYIQLlmeIHzRybQskk9nygiuolRsTNwI+F5ldU5ODgkJCWRmZlodigCCgoKoX78+gYGBpXZOSRAu2OaglgpqUQoe/fVRwPfqIBISEqhcuTLh4eEomU3RUlprjh8/TkJCAhEREaV2Xm82c/0KiALClFIJwDNa60+9dT5PiksyWjA1qSUlCOF970W/Z3UILmVmZkpy8BFKKapXr86xY8dK9bzebMV0q7eO7W1xh1OpWy2YikGlV5QTZVebmr7bTUiSg++w4m8hYzG5EJeUKh3kRKlZc3ANaw6usToMIZxIgijkdEYOSSkZkiBEqXly+ZM8ufxJq8MQLqxYsYKBAwc6rY+NjWXRokXndcwXXzw7vHt8fDxt2vhuCVISRCH7bBXUUv8gSsdHAz/io4EfWR3GRSs3N7fUz1lcgigpHvsE4eukFVMhewsShJQgRClpHtbc6hDcExXlvO7mm2HcOEhPh+ho5+2jRxuP5GQYNsxx24oVJZ5y6tSpzJkzhwYNGhAWFkbnzp2ZOHEiUVFR9OjRg9WrVzN48GCGDRvGnXfeybFjx6hRowaff/45DRs2ZPTo0QwcOJBh5rkrVapEWloaK1as4NlnnyUsLIwdO3bQuXNnZs+ejVKKJUuWMGHCBMLCwujUqZNTTNnZ2Tz99NNkZGSwatUqJk2axK5duzh06BDx8fGEhYVx7bXXEhMTw3vvGQ0QBg4cyMSJE1myZAkZGRl06NCB1q1b89///pe8vDzGjBnDmjVrqFevHvPnz6dChQolXpvSICWIQuIOn6JmSAVCgstZHYooI1bGr2Rl/Eqrw/A5MTExfP/992zZsoV58+YRExPjsD0lJYWVK1fyyCOPMH78eG677Ta2bdvGiBEjeOCBB0o8/pYtW3j77bfZuXMn+/fvZ/Xq1WRmZjJmzBgWLFjAn3/+SVJSktPrypUrx/PPP8/w4cOJjY1l+PDhAGzatIn58+fz5ZdfFnnOl19+mQoVKhAbG8ucOXMA2Lt3L/fddx9//fUXoaGhfP/99+dymbxKShCFGBXUcntJlJ5nVjwD+F4/CCfF/eIPDi5+e1iYWyUGe6tWrWLIkCG2X9ODBg1y2F7wxQywdu1a5s2bB8CoUaN47LHHSjz+5ZdfTv369QHo0KED8fHxVKpUiYiICJo1awbAyJEjmT59ulvxDh48+Lx++UdERNChg9GbvnPnzsTHx5/zMbxFEoSdM5k5JB4/w9Vt6lkdiihDPhvymdUh+CStix95p2LFikVuK2gSGhAQQH5+vu142dnZtn3Kly9vW/b397fVHZxvc1L7eOzPCxTbG71wHBkZGed1fm+QW0x2dh9KQQMt6odaHYooQxpXbUzjqo2tDsPn9OrViwULFpCZmUlaWhoLFy4sct8ePXrw9ddfAzBnzhx69eoFQHh4OJs2bQJg/vz55OTkFHvOFi1acODAAfbt2wfAV1995XK/ypUrc/r06SKPEx4eTmxsLPn5+Rw8eJANGzbYtgUGBpYYh6+QBGFnV0IKCmhRTxKEKD3L9i9j2f5lVofhc7p06cLgwYNp3749Q4cOJTIykpAQ141Hpk2bxueff067du2YNWsW77zzDgBjxoxh5cqVXH755axfv77YUgcY4x1Nnz6dAQMG0KtXLxo1auRyv969e7Nz5046dOjAN99847S9Z8+eRERE0LZtWyZOnOhQ2T127FjatWvHiBEj3L0UllElFeNKU2RkpC5cEVWaJn+5geRTmXz0nysti0GUPVEzogDfq4PYtWsXLVu2tDSGtLQ0KlWqRHp6OldeeSXTp0932bKorHD1N1FKbdJaR3rjfFIHYcrXmt2JJ+nVso7VoYgyZtYNs6wOwWeNHTuWnTt3kpmZye23316mk4MVJEGYEpLTSMvMpVX9qlaHIsqYBiENrA7BZxXXZFR4n9RBmHYlpgDQUuofRClbEreEJXFLrA5DCCdSgjDtSjhJpaAA6odVsjoUUca8vOplAPo37W9xJEI4kgRh2vHvCVrUq4qfDG8sStnXw762OgQhXJIEARxNzeDg8TNc16mh1aGIMqh2pdpWhyCES1IHAcTsM2ZpimxSw+JIRFm0YM8CFuxZYHUYPumtt96idevWtGnThltvvdXWI/nAgQN07dqVZs2aMXz4cFsP6XfffZc2bdoQHR1tW7dq1SoefvjhUo27qGHC3eUrw4BLgsBIEDWqBNFQ6h+EBd5Y+wZvrH3D6jB8TmJiItOmTSMmJoYdO3aQl5dn6y39+OOP89BDD7F3716qVq3Kp58asxl/8sknbNu2jY4dO7J06VK01kydOpUpU6aUeD6ttcPwGEISBLl5+Ww5kExkkxoyvaKwxNyb5zL35rlWh1GiqBlRzIidAUBOXg5RM6KYvW02AOk56UTNiOKbHUav4tTMVKJmRDFvlzGAXnJ6MlEzomwlpaQ051FSXcnNzSUjI4Pc3FzS09OpW7cuWmt+++032xDet99+Oz/++KPtNTk5OaSnpxMYGMisWbOIjo6malXXzdfj4+Np2bIl48aNo1OnThw8eJBffvmF7t2706lTJ2666SbS0tIAeP755+nSpQtt2rRh7NixtrGi4uLi6NOnD+3bt6dTp062YTrS0tIYNmwYLVq0YMSIEbb9N23axFVXXUXnzp3p168fhw8ftq1v37493bt353//+59b18fbynyC2J2YQnpWLp3l9pKwSFhwGGHBYVaH4XPq1avHxIkTadiwIXXq1CEkJIRrr72W48ePExoaSkCAUYVav359EhMTAZg4cSLdunXj2LFj9OzZk5kzZzJu3Lhiz7Nnzx5uu+02tmzZQsWKFXnhhRdYtmwZmzdvJjIykjfffBOA8ePHs3HjRnbs2EFGRgY///wzACNGjOC+++5j69atrFmzhjp1jM62roYTz8nJ4f7772fu3Lls2rSJO++8k8mTJwNwxx13MG3aNNauXeuV63k+ynwldcy+Y/gpRccI+Q8qrFHwK3toy6EWR1I8+6FAAv0DHZ4HBwY7PA8JCnF4HhYc5vDcnYr5kydPMn/+fA4cOEBoaCg33XQTs2fPpl+/fk77FpT+R40axahRowB47rnneOCBB1i8eDFffPEFDRo04I033sDPz/F3caNGjejWrRsA69atY+fOnfTs2RMwJgfq3r07AL///juvvvoq6enpnDhxgtatWxMVFUViYiI33HADYIzlVMDVcOKhoaHs2LGDvn37ApCXl0edOnVITU0lJSWFq666yvY+Fi9eXOI18rYynSC01vyx8zBtG1WjUlCg1eGIMmra+mmA7yeI0rZs2TIiIiKoUcMo3Q8dOpQ1a9YwYsQIUlJSyM3NJSAggISEBOrWrevw2kOHDrFx40aeeeYZLr/8ctauXcvkyZNZvny57cu5gP0Aflpr+vbt6zSKa2ZmJuPGjSMmJoYGDRrw7LPPkpmZWeyQ5K6GE9da07p1a6dSQkpKik/e4i7Tt5j2HEol8cQZrmkr8z8I68y/ZT7zb5lvdRg+p2HDhqxbt4709HS01ixfvpyWLVuilKJ3797MnWvU28ycOZMhQ4Y4vHbKlClMnToVgIyMDJRS+Pn5kZ6eXuw5u3XrxurVq4mLiwMgPT2dv//+29Z6KiwsjLS0NNu5q1SpQv369W11IFlZWcWeo3nz5hw7dsyWIHJycmwzyYWEhLBq1SoA22xzVivTCeK37YkE+vvRq4W0QxfWCQkKISRI5kAvrGvXrgwbNoxOnTrRtm1b8vPzGTt2LACvvPIKb775Jk2bNuX48ePcddddttdt2bIFgI4dOwJw11130bZtWzZv3kz//sX3Vq9RowYzZszg1ltvpV27dnTr1o3du3cTGhrKmDFjaNu2Lddffz1dunSxvWbWrFlMmzaNdu3a0aNHD5fTlBYoV64cc+fO5fHHH6d9+/Z06NCBNWvWAPD5559z33330b17d5+Zk7rMDvd9JjOHkdN+4/KmNZk0tGOpnFMIVwpa/gxvM7yEPUuXLwz3LRzJcN+lZNHmf0nPymVYd5nJS1jrg5gPAN9LEEKUyQRxJjOHuev20yGiOs3qSNFeWGvRiEVWhyCES2UyQXy5Ko7UM9ncdXULq0MRguDAYKtDKJLW2idb15RFVlQHlLlK6r8OnmDeugP069CAy+rK3A/CerO3zbb1SPYlQUFBHD9+3JIvJuFIa83x48cd+lmUhjJVgvj7UApPf72RWqEVGNtXKt+Eb/hk8ycAjGw30uJIHNWvX5+EhASOHTtmdSgCI2EXdLwrLWUmQazZncTrP22lUoVAXhnZlYrSMU74iF9H/Wp1CC4FBgYSERFhdRjCQl5NEEqp/sA7gD/widb6ZW+er7CcvHxi4o6xICaeTfuTaVq7Cs/cHEnNEN9oYywEGMNWCOGLvJYglFL+wP+AvkACsFEp9ZPWeue5HCdfa7TW5OVrtDae52tNbp4mMzuXjOw8MnPyyMzOJTU9m+TTmRxNzWD/kVPsPZxKRnYeoRXLcfc1Lbi+awSB/mWu2kX4uIIRUkd3GG1pHEIU5s0SxOVAnNZ6P4BS6mtgCFBkgohLOsXAFxcb47JrTf551o2VD/QnomZl+ravz+VNa9IxIowASQzCR0mCEL7KmwmiHnDQ7nkC0LXwTkqpscBY82nWwsnRO7wYkyeEAclWB+EGidOzvB6nusMjzUnlenrWxRBnc28d2JsJwtWn3alMoLWeDkwHUErFeKvLuKdcDDGCxOlpEqdnSZyeo5Ty2vhE3rzvkgA0sHteHzjkxfMJIYTwIG8miI1AM6VUhFKqHHAL8JMXzyeEEMKDvHaLSWudq5QaDyzFaOb6mdb6rxJeNt1b8XjQxRAjSJyeJnF6lsTpOV6L0aeG+xZCCOE7pO2nEEIIlyRBCCGEcMknEoRSqr9Sao9SKk4p9YQF52+glPpdKbVLKfWXUupBc/2zSqlEpVSs+Yi2e80kM949Sql+pfVelFLxSqntZjwx5rpqSqlflVJ7zX+rmuuVUmqaGcs2pVQnu+Pcbu6/Vyl1uwfja253vWKVUqeUUhN84VoqpT5TSh1VSu2wW+exa6eU6mz+beLM155Xx4Yi4nxNKbXbjOUHpVSouT5cKZVhd10/LCmeot6zh+L02N9ZGQ1c1ptxfqOMxi6eivMbuxjjlVKx5npLrqcq+jvI2s+nNoeysOqBUYG9D2gMlAO2Aq1KOYY6QCdzuTLwN9AKeBaY6GL/Vmac5YEIM37/0ngvQDwQVmjdq8AT5vITwCvmcjSwGKNPSjdgvbm+GrDf/LequVzVS3/bJKCRL1xL4EqgE7DDG9cO2AB0N1+zGLjOg3FeCwSYy6/YxRluv1+h47iMp6j37KE4PfZ3Br4FbjGXPwTu9VSchba/ATxt5fWk6O8gSz+fvlCCsA3JobXOBgqG5Cg1WuvDWuvN5vJpYBdGT/CiDAG+1lpnaa0PAHEY78Oq9zIEmGkuzwSut1v/hTasA0KVUnWAfsCvWusTWuuTwK9A8bO5n59rgH1a639KiL1UrqXW+g/ghIvzX/C1M7dV0Vqv1cb/xi/sjnXBcWqtf9Fa55pP12H0KypSCfEU9Z4vOM5inNPf2fx1ezUw15txmue5GfiquGN4+3oW8x1k6efTFxKEqyE5ivty9iqlVDjQEVhvrhpvFuE+sys6FhVzabwXDfyilNqkjGFKAGpprQ+D8UEDavpAnGD0fbH/j+dr1xI8d+3qmcvejhfgToxfgAUilFJblFIrlVJXmOuKi6eo9+wpnvg7VwdS7JKit67nFcARrfVeu3WWXs9C30GWfj59IUG4NSRHaVBKVQK+ByZorU8BHwBNgA7AYYyiKBQdc2m8l55a607AdcB9Sqkri9nXsjjN+8WDge/MVb54LYtzrnGVSrxKqclALjDHXHUYaKi17gg8DHyplKpSWvG44Km/c2nFfyuOP2IsvZ4uvoOK3LWIeDx6PX0hQfjEkBxKqUCMP8wcrfU8AK31Ea11ntY6H/gYozgMRcfs9feitT5k/nsU+MGM6YhZhCwoCh+1Ok6MBLZZa33EjNfnrqXJU9cuAcfbPh6P16xwHAiMMG8TYN6yOW4ub8K4n39ZCfEU9Z4vmAf/zskYt00CCq33GPPYQ4Fv7OK37Hq6+g4q5til8/k818oUTz8wenPvx6i4Kqikal3KMSiMe3JvF1pfx275IYx7qACtcaxw249R2ebV9wJUBCrbLa/BqDt4DceKrFfN5QE4VmRt0Gcrsg5gVGJVNZerefiafg3c4WvXkkKVkJ68dhjDy3TjbCVgtAfj7I8xVH6NQvvVAPzN5cZAYknxFPWePRSnx/7OGKVP+0rqcZ6K0+6arvSF60nR30GWfj499oVwIQ+MGvm/MbL1ZAvO3wujuLUNiDUf0cAsYLu5/qdCH/7JZrx7sGsN4M33Yn5gt5qPvwqOj3G/djmw1/y34AOhMCZt2me+j0i7Y92JUVEYh90XuYfiDAaOAyF26yy/lhi3Eg4DORi/qO7y5LUDIoEd5mvewxypwENxxmHcWy74fH5o7nuj+VnYCmwGBpUUT1Hv2UNxeuzvbH7eN5jv/TugvKfiNNfPAP5TaF9LridFfwdZ+vmUoTaEEEK45At1EEIIIXyQJAghhBAuSYIQQgjhkiQIIYQQLkmCEEII4ZIkCHFRUkqtUEp5fTJ5pdQD5gibcwqtj1RKTTOXo5RSPTx4znCl1P+5OpcQpclrU44K4auUUgH67Bg/JRmH0Wb/gP1KrXUMEGM+jQLSMDoueiKGcOD/gC9dnEuIUiMlCOE15i/hXUqpj80x7n9RSlUwt9lKAEqpMKVUvLk8Win1o1JqgVLqgFJqvFLqYXPwtHVKqWp2pxiplFqjlNqhlLrcfH1Fc5C4jeZrhtgd9zul1ALgFxexPmweZ4dSaoK57kOMzlo/KaUeKrR/lFLqZ3Ngtf8ADylj/oArlFI1lFLfmzFsVEr1NF/zrFJqulLqF+AL8/r8qZTabD4KSiEvA1eYx3uo4FzmMaqZ12ebeT3a2R37M/O67ldKPWB3PRYqpbaa7234hf1VRZlyob1T5SGPoh4Yv4RzgQ7m82+BkebyCszen0AYEG8uj8boAVoZY9iDVMzersBbGIOYFbz+Y3P5SsxhFIAX7c4RitFDt6J53ARc9HIFOmP0Rq0IVMLoSdvR3BZPofk3zPVRwM/m8rPYzYGA8cu/l7ncENhlt98moIL5PBgIMpebATGFj+3iXO8Cz5jLVwOxdsdegzGURRhGT/ZAjJ7BH9sdK6Twe5GHPIp6yC0m4W0HtNax5vImjKRRkt+1MSb+aaVUKrDAXL8daGe331dgjPevlKqijFnWrgUGK6UmmvsEYXxJgzlOvovz9QJ+0FqfAVBKzcMYBnqLO2/QhT5AK3V2wq4qSqnK5vJPWusMczkQeE8p1QHIwxgUriS9ML700Vr/ppSqrpQKMbct1FpnAVlKqaNALYxr9rpS6hWMJPPneb4nUQZJghDelmW3nAdUMJdzOXuLM6iY1+TbPc/H8TNbeJyYgmGNb9Ra77HfoJTqCpwpIsbzmhq0GH5Ad7tEUBADhWJ4CDgCtDdfk+nGsYsbtrnwtQ7QWv+tlOqMMa7PS0qpX7TWz7v1LkSZJ3UQwirxGLd2AIad5zGGAyilegGpWutUYClwv1K2+YI7unGcP4DrlVLBSqmKwA3AufzSPo1xS6zAL8D4gidmCcGVEOCwNobGHoUxuqmr4xWOdYR53CggWRczb4BSqi6QrrWeDbyOMfWmEG6RBCGs8jpwr1JqDcY98/Nx0nz9hxgjiQJMxbh1s00Zk9RPLekg2pjqcQbGyKHrgU+01udye2kBcENBJTXwABBpViTvxKjEduV94Hal1DqM20sFpYttQK5ZsfxQodc8W3BsjMrs20uIrS2wQSkVizGa6gvn8L5EGSejuQohhHBJShBCCCFckgQhhBDCJUkQQgghXJIEIYQQwiVJEEIIIVySBCGEEMIlSRBCCCFc+n83txowrIKr6gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "for i in range(rep):\n",
    "    plt.plot(mi_list[i,:],color='steelblue')\n",
    "    for t in range(mi_list[i].shape[0]):\n",
    "        if (mi_list[0,t]>.8*mi):\n",
    "            plt.axvline(t,label='80% reached',linestyle=':',color='green')\n",
    "            break\n",
    "plt.xlim((0,mi_list[0].shape[0]))\n",
    "plt.ylim((0,mi*1.1))\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig(fig_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
